EDA

Set-up

library(tidyverse)
library(janitor)
library(readxl)
library(writexl)
library(patchwork)
library(glmmTMB)
library(car)
library(emmeans)
library(nlme)
library(lmerTest)
library(DHARMa)
library(brms)

theme_set(theme_minimal())

Data Read-in

Convert 3 .csv to 1 Excel File (DO NOT NEED TO RUN)

df_combined <- read_csv(here::here("data", "Combined_Master.csv")) %>% clean_names()
df_year_1 <- read_csv(here::here("data", "Year1_Master.csv")) %>% clean_names()
df_year_2 <- read_csv(here::here("data", "Year2_Master.csv")) %>% clean_names()

# Change date format
df_combined <- df_combined %>%
  mutate(date = as.Date(date, format="%m_%d_%y")) %>%
  select(-starts_with("x"))

df_year_1 <- df_year_1 %>%
  mutate(date = as.Date(date, format="%m_%d_%y")) %>%
  select(-starts_with("x"))

df_year_2 <- df_year_2 %>%
  mutate(date = as.Date(date, format="%m_%d_%y")) %>%
  select(-starts_with("x"))

# Save as single Excel File
out <- list("combined" = df_combined, "year_1" = df_year_1, "year_2" = df_year_2)
write_xlsx(out, "data/Master.xlsx")

Read-in Data from Excel File

df_combined <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "combined")
df_year_1 <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "year_1")
df_year_2 <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "year_2")

df_combined <- df_combined %>%
  mutate(
    sample_id = as.factor(sample_id),
    site = as.factor(site),
    block = as.factor(block),
    study_year = as.factor(study_year)
  )

df_year_1 <- df_year_1 %>%
  mutate(
    sample_id = as.factor(sample_id),
    site = as.factor(site),
    block = as.factor(block),
    study_year = as.factor(study_year)
  )

df_year_2 <- df_year_2 %>%
  mutate(
    sample_id = as.factor(sample_id),
    site = as.factor(site),
    block = as.factor(block),
    study_year = as.factor(study_year)
  )

EDA

Create Datasets where LOD is 0 or NA

response_vars <- c("nh4", "no3", "mg", "p", "ec", "p_h")

df_combined_LOD_0 <- df_combined %>%
    mutate(across(
    all_of(response_vars), 
    ~ {
      # Replace "LOD" with 0
      x <- ifelse(.x == "LOD", 0, .x)
      as.numeric(x)
    }
  ))

df_combined_LOD_NA <- df_combined %>%
    mutate(across(
    all_of(response_vars), 
    ~ {
      # Replace "LOD" with NA
      x <- ifelse(.x == "LOD", NA, .x)
      as.numeric(x)
    }
  ))

Count of LOD

df_combined %>%
  group_by(treatment) %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x == "LOD", na.rm = T)
      }
    ))

# A tibble: 6 × 7
  treatment   nh4   no3      mg     p    ec   p_h
  <chr>     <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>
1 MC0       0.774 0.708 0       0.317     0     0
2 MC1       0.843 0.746 0       0.405     0     0
3 MC2       0.848 0.741 0       0.431     0     0
4 PC0       0.850 0.685 0       0.248     0     0
5 PC1       0.776 0.724 0.00971 0.186     0     0
6 PC2       0.836 0.694 0.00962 0.176     0     0

Count of Zeros

prop_zero_0 <- df_combined_LOD_0 %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x == 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = everything(), names_to = "response", values_to = "prop_zero")

prop_zero_NA <- df_combined_LOD_NA %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x == 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = everything(), names_to = "response", values_to = "prop_zero")

prop_zero_0 %>%
  left_join(prop_zero_NA, by = "response", suffix = c("_0", "_NA")) %>%
  arrange(-prop_zero_0)

# A tibble: 6 × 3
  response prop_zero_0 prop_zero_NA
  <chr>          <dbl>        <dbl>
1 nh4          0.822              0
2 no3          0.716              0
3 p            0.297              0
4 mg           0.00316            0
5 ec           0                  0
6 p_h          0                  0

prop_zero_0 <- df_combined_LOD_0 %>%
  group_by(study_year) %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x == 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = -study_year, names_to = "response", values_to = "prop_zero")

prop_zero_NA <- df_combined_LOD_NA %>%
  group_by(study_year) %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x == 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = -study_year, names_to = "response", values_to = "prop_zero")

prop_zero_0 %>%
  left_join(prop_zero_NA, by = c("response", "study_year"), suffix = c("_0", "_NA")) %>%
  arrange(study_year, -prop_zero_0)

# A tibble: 12 × 4
   study_year response prop_zero_0 prop_zero_NA
   <fct>      <chr>          <dbl>        <dbl>
 1 1          nh4          0.842              0
 2 1          no3          0.709              0
 3 1          p            0.192              0
 4 1          mg           0.0025             0
 5 1          ec           0                  0
 6 1          p_h          0                  0
 7 2          nh4          0.789              0
 8 2          no3          0.728              0
 9 2          p            0.484              0
10 2          mg           0.00429            0
11 2          ec           0                  0
12 2          p_h          0                  0

Negative Values

prop_neg_0 <- df_combined_LOD_0 %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x < 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = everything(), names_to = "response", values_to = "prop_negative")

prop_neg_NA <- df_combined_LOD_NA %>%
  summarise(across(
    all_of(response_vars), 
    ~ {
      x <- mean(.x < 0, na.rm = T)
      }
    )) %>%
  pivot_longer(cols = everything(), names_to = "response", values_to = "prop_negative")

prop_neg_0 %>%
  left_join(prop_neg_NA, by = "response", suffix = c("_0", "_NA")) %>%
  arrange(-prop_negative_0)

# A tibble: 6 × 3
  response prop_negative_0 prop_negative_NA
  <chr>              <dbl>            <dbl>
1 nh4                    0                0
2 no3                    0                0
3 mg                     0                0
4 p                      0                0
5 ec                     0                0
6 p_h                    0                0

Concentration by Treatment Box Plots

capwords <- function(s, strict = FALSE) {
    cap <- function(s) paste(toupper(substring(s, 1, 1)),
                  {s <- substring(s, 2); if(strict) tolower(s) else s},
                             sep = "", collapse = " " )
    sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))
}

plot_concentration <- function(r, data) {
  data %>% 
    filter(.data[[r]] > 0) %>%
    ggplot(aes(x = treatment, y = .data[[r]], fill = treatment)) +
      geom_boxplot(alpha = 0.7) +
      geom_point(alpha = 0.3) +
      facet_wrap(vars(study_year), scales = "free") +
      labs(title = paste0(capwords(r), " Concentration (Non-Zero) by Treatment"), x = "Treatment", y = capwords(r)) +
      theme_bw()
}

plots_0 <- map(response_vars, plot_concentration, data = df_combined_LOD_0)
wrap_plots(plots_0, ncol = 2, guides = "collect")

# It will be the same since we are plotting non-zeros
# plots_NA <- map(response_vars, plot_concentration, data = df_combined_LOD_NA)
# wrap_plots(plots_NA, ncol = 2, guides = "collect")

Plot Response Variables Over Time

capwords <- function(s, strict = FALSE) {
    cap <- function(s) paste(toupper(substring(s, 1, 1)),
                  {s <- substring(s, 2); if(strict) tolower(s) else s},
                             sep = "", collapse = " " )
    sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))
}

plot_response <- function(r, data, by) {
  data %>%
    ggplot(aes(x = date, y = .data[[r]], color = .data[[by]])) +
    facet_wrap(vars(study_year), scales = "free") +
    geom_point() +
    geom_line() +
    theme_bw(base_size = 12) +
    theme(legend.position = "none") +
    labs(
      title = paste0(capwords(r), " Over Time"),
      x = "Date"
    )
}
  
plots_0 <- map(response_vars, plot_response, data = df_combined_LOD_0, by = "sample_id")
wrap_plots(plots_0, ncol = 2)

plots_NA <- map(response_vars, plot_response, data = df_combined_LOD_NA, by = "sample_id")
wrap_plots(plots_NA, ncol = 2)

Models

Responses with no Zero

Data

df_model <- df_combined_LOD_0 %>%
  filter(study_year == 1) %>% # Can change this to only look at year 2, or comment out to look at both years
  mutate(
    date = factor(date),
    treatment = factor(treatment),
    site = factor(site),
    block = factor(block),
    
    # Create combined grouping factors
    site_block_id = interaction(site, block, drop = TRUE),
    trt_site_block_id = interaction(treatment, site, block, drop = TRUE)
  )

EC

# fit_lme_ec <- lme(
#   fixed = ec ~ date * treatment,
#   random = list(
#     site_block_id = ~ 1, # Random intercept for site:block combinations
#     trt_site_block_id = ~ 1 # Random intercept for treatment:site:block combinations
#   ),
#   data = df_model,
#   na.action = na.exclude
# )
# 
# Anova(fit_lme_ec)
# 
# plot(fitted(fit_lme_ec), residuals(fit_lme_ec),
#      xlab = "Fitted Values",
#      ylab = "Residuals",
#      main = "Residuals vs. Fitted Values")
# abline(h = 0, col = "red", lty = 2)
# 
# emmip(fit_lme_ec, treatment ~ date, type = "response", CIs = T) + theme_bw()

fit_glmer_gamma_ec <- glmer(
  ec ~ date * treatment +
    (1 | site_block_id) +
    (1 | trt_site_block_id),
  data = df_model,
  family = Gamma(link = "log"),
  control = glmerControl(optimizer = "bobyqa"), # default optimizer did not converge
  na.action = na.exclude
)

Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
length(par)^2 is not recommended.

summary(fit_glmer_gamma_ec)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: ec ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
   Data: df_model
Control: glmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
   147.1    375.2    -16.6     33.1      347 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9861 -0.5108 -0.0331  0.4374  4.3962 

Random effects:
 Groups            Name        Variance Std.Dev.
 trt_site_block_id (Intercept) 0.03329  0.1825  
 site_block_id     (Intercept) 0.02424  0.1557  
 Residual                      0.06125  0.2475  
Number of obs: 404, groups:  trt_site_block_id, 48; site_block_id, 8

Fixed effects:
                             Estimate Std. Error t value Pr(>|z|)    
(Intercept)                  0.276677   0.235992   1.172  0.24104    
date2024-01-12              -0.029280   0.125662  -0.233  0.81576    
date2024-01-23              -0.155530   0.121384  -1.281  0.20009    
date2024-02-07              -0.123118   0.121340  -1.015  0.31027    
date2024-02-21              -0.077465   0.125634  -0.617  0.53750    
date2024-03-06              -0.109187   0.121636  -0.898  0.36937    
date2024-03-22               0.018036   0.125707   0.143  0.88592    
date2024-04-04              -0.104127   0.131053  -0.795  0.42688    
date2024-04-16              -0.268501   0.124017  -2.165  0.03038 *  
treatmentMC1                 0.351553   0.218209   1.611  0.10716    
treatmentMC2                 0.345026   0.221048   1.561  0.11856    
treatmentPC0                -0.374132   0.218786  -1.710  0.08726 .  
treatmentPC1                -0.027530   0.217548  -0.127  0.89930    
treatmentPC2                -0.246458   0.229081  -1.076  0.28199    
date2024-01-12:treatmentMC1 -0.084005   0.168474  -0.499  0.61804    
date2024-01-23:treatmentMC1  0.005998   0.163360   0.037  0.97071    
date2024-02-07:treatmentMC1 -0.238434   0.163187  -1.461  0.14399    
date2024-02-21:treatmentMC1 -0.364844   0.166539  -2.191  0.02847 *  
date2024-03-06:treatmentMC1 -0.469824   0.163802  -2.868  0.00413 ** 
date2024-03-22:treatmentMC1 -0.450368   0.166540  -2.704  0.00685 ** 
date2024-04-04:treatmentMC1 -0.471514   0.173754  -2.714  0.00665 ** 
date2024-04-16:treatmentMC1 -0.609815   0.165547  -3.684  0.00023 ***
date2024-01-12:treatmentMC2 -0.015509   0.169959  -0.091  0.92729    
date2024-01-23:treatmentMC2  0.173799   0.166840   1.042  0.29755    
date2024-02-07:treatmentMC2 -0.080414   0.167165  -0.481  0.63048    
date2024-02-21:treatmentMC2 -0.041221   0.170343  -0.242  0.80879    
date2024-03-06:treatmentMC2  0.020059   0.167565   0.120  0.90471    
date2024-03-22:treatmentMC2  0.005482   0.170444   0.032  0.97434    
date2024-04-04:treatmentMC2  0.006028   0.174387   0.035  0.97243    
date2024-04-16:treatmentMC2 -0.147860   0.172094  -0.859  0.39024    
date2024-01-12:treatmentPC0 -0.063463   0.166513  -0.381  0.70311    
date2024-01-23:treatmentPC0 -0.107943   0.163261  -0.661  0.50851    
date2024-02-07:treatmentPC0  0.096256   0.163313   0.589  0.55559    
date2024-02-21:treatmentPC0  0.092224   0.166460   0.554  0.57956    
date2024-03-06:treatmentPC0 -0.103295   0.163453  -0.632  0.52742    
date2024-03-22:treatmentPC0 -0.034863   0.166469  -0.209  0.83411    
date2024-04-04:treatmentPC0 -0.161729   0.170559  -0.948  0.34301    
date2024-04-16:treatmentPC0 -0.117894   0.165458  -0.713  0.47614    
date2024-01-12:treatmentPC1 -0.214503   0.168533  -1.273  0.20310    
date2024-01-23:treatmentPC1 -0.182493   0.165601  -1.102  0.27046    
date2024-02-07:treatmentPC1 -0.298994   0.163356  -1.830  0.06720 .  
date2024-02-21:treatmentPC1 -0.274619   0.166537  -1.649  0.09915 .  
date2024-03-06:treatmentPC1 -0.390718   0.163766  -2.386  0.01704 *  
date2024-03-22:treatmentPC1 -0.327529   0.166667  -1.965  0.04939 *  
date2024-04-04:treatmentPC1 -0.438312   0.170963  -2.564  0.01035 *  
date2024-04-16:treatmentPC1 -0.545116   0.170793  -3.192  0.00141 ** 
date2024-01-12:treatmentPC2  0.051349   0.187448   0.274  0.78413    
date2024-01-23:treatmentPC2 -0.108320   0.178830  -0.606  0.54470    
date2024-02-07:treatmentPC2 -0.052320   0.179043  -0.292  0.77012    
date2024-02-21:treatmentPC2 -0.035076   0.181874  -0.193  0.84707    
date2024-03-06:treatmentPC2 -0.190791   0.179275  -1.064  0.28722    
date2024-03-22:treatmentPC2 -0.059831   0.181894  -0.329  0.74221    
date2024-04-04:treatmentPC2 -0.187168   0.187910  -0.996  0.31922    
date2024-04-16:treatmentPC2 -0.265369   0.180586  -1.469  0.14170    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 54 > 12.
Use print(value, correlation=TRUE)  or
    vcov(value)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
maxfun < 10 * length(par)^2 is not recommended.

plot(fitted(fit_glmer_gamma_ec), residuals(fit_glmer_gamma_ec),
     xlab = "Fitted Values",
     ylab = "Residuals",
     main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)

emmip(fit_glmer_gamma_ec, treatment ~ date, type = "response", CIs = T) + theme_bw()

Ph

fit_lme_ph <- lme(
  fixed = p_h ~ as.factor(date) * treatment,
  random = list(
    site_block_id = ~ 1, # Random intercept for site:block combinations
    trt_site_block_id = ~ 1 # Random intercept for treatment:site:block combinations
  ),
  data = df_model,
  na.action = na.exclude
)

Anova(fit_lme_ph)

Analysis of Deviance Table (Type II tests)

Response: p_h
                           Chisq Df Pr(>Chisq)    
as.factor(date)           59.922  8  4.829e-10 ***
treatment                 24.529  5  0.0001717 ***
as.factor(date):treatment 20.744 40  0.9949031    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

S(fit_lme_ph)

Linear mixed model fit by REML,  Data: df_model

Fixed Effects:
 Formula: p_h ~ as.factor(date) * treatment 

                                         Estimate  Std.Error  df t value
(Intercept)                             8.0741630  0.1409290 302  57.292
as.factor(date)2024-01-12               0.1075456  0.1441116 302   0.746
as.factor(date)2024-01-23               0.0358370  0.1394633 302   0.257
as.factor(date)2024-02-07               0.0995870  0.1394633 302   0.714
as.factor(date)2024-02-21              -0.1111751  0.1440079 302  -0.772
as.factor(date)2024-03-06              -0.0391630  0.1394633 302  -0.281
as.factor(date)2024-03-22              -0.1610258  0.1441116 302  -1.117
as.factor(date)2024-04-04              -0.0853503  0.1499612 302  -0.569
as.factor(date)2024-04-16              -0.1474028  0.1427589 302  -1.033
treatmentMC1                           -0.4743938  0.1612606  35  -2.942
treatmentMC2                           -0.4064689  0.1612408  35  -2.521
treatmentPC0                           -0.1746424  0.1564647  35  -1.116
treatmentPC1                           -0.4088924  0.1564907  35  -2.613
treatmentPC2                           -0.3024218  0.1771349  35  -1.707
as.factor(date)2024-01-12:treatmentMC1  0.2313875  0.1975466 302   1.171
as.factor(date)2024-01-23:treatmentMC1  0.4002817  0.1981007 302   2.021
as.factor(date)2024-02-07:treatmentMC1  0.2643938  0.1917050 302   1.379
as.factor(date)2024-02-21:treatmentMC1  0.2601559  0.1950219 302   1.334
as.factor(date)2024-03-06:treatmentMC1  0.2893938  0.1917050 302   1.510
as.factor(date)2024-03-22:treatmentMC1  0.3512565  0.1950973 302   1.800
as.factor(date)2024-04-04:treatmentMC1  0.2825575  0.2025085 302   1.395
as.factor(date)2024-04-16:treatmentMC1  0.2526336  0.1941574 302   1.301
as.factor(date)2024-01-12:treatmentMC2 -0.0064897  0.1950813 302  -0.033
as.factor(date)2024-01-23:treatmentMC2  0.0794538  0.1980666 302   0.401
as.factor(date)2024-02-07:treatmentMC2  0.1427189  0.1916883 302   0.745
as.factor(date)2024-02-21:treatmentMC2  0.0697310  0.1950072 302   0.358
as.factor(date)2024-03-06:treatmentMC2  0.0377189  0.1916883 302   0.197
as.factor(date)2024-03-22:treatmentMC2  0.0395817  0.1950813 302   0.203
as.factor(date)2024-04-04:treatmentMC2  0.1551562  0.1994248 302   0.778
as.factor(date)2024-04-16:treatmentMC2 -0.0322732  0.1971682 302  -0.164
as.factor(date)2024-01-12:treatmentPC0 -0.1445663  0.1911626 302  -0.756
as.factor(date)2024-01-23:treatmentPC0 -0.0991076  0.1876887 302  -0.528
as.factor(date)2024-02-07:treatmentPC0 -0.1791076  0.1876887 302  -0.954
as.factor(date)2024-02-21:treatmentPC0 -0.0870955  0.1910851 302  -0.456
as.factor(date)2024-03-06:treatmentPC0 -0.1316076  0.1876887 302  -0.701
as.factor(date)2024-03-22:treatmentPC0 -0.0959949  0.1911626 302  -0.502
as.factor(date)2024-04-04:treatmentPC0 -0.0754203  0.1956039 302  -0.386
as.factor(date)2024-04-16:treatmentPC0 -0.2146178  0.1901452 302  -1.129
as.factor(date)2024-01-12:treatmentPC1  0.0724544  0.1936173 302   0.374
as.factor(date)2024-01-23:treatmentPC1  0.2155916  0.1901828 302   1.134
as.factor(date)2024-02-07:treatmentPC1  0.2263924  0.1877103 302   1.206
as.factor(date)2024-02-21:treatmentPC1  0.2259045  0.1911024 302   1.182
as.factor(date)2024-03-06:treatmentPC1  0.1713924  0.1877103 302   0.913
as.factor(date)2024-03-22:treatmentPC1  0.3495051  0.1911795 302   1.828
as.factor(date)2024-04-04:treatmentPC1  0.2163297  0.1956153 302   1.106
as.factor(date)2024-04-16:treatmentPC1  0.1436748  0.1964970 302   0.731
as.factor(date)2024-01-12:treatmentPC2 -0.0387264  0.2149142 302  -0.180
as.factor(date)2024-01-23:treatmentPC2  0.1974218  0.2052379 302   0.962
as.factor(date)2024-02-07:treatmentPC2  0.0036718  0.2052379 302   0.018
as.factor(date)2024-02-21:treatmentPC2  0.1006839  0.2083860 302   0.483
as.factor(date)2024-03-06:treatmentPC2  0.0672609  0.2074953 302   0.324
as.factor(date)2024-03-22:treatmentPC2  0.1380346  0.2084665 302   0.662
as.factor(date)2024-04-04:treatmentPC2  0.1142633  0.2147977 302   0.532
as.factor(date)2024-04-16:treatmentPC2  0.0006617  0.2074494 302   0.003
                                       Pr(>|t|)    
(Intercept)                             < 2e-16 ***
as.factor(date)2024-01-12               0.45609    
as.factor(date)2024-01-23               0.79738    
as.factor(date)2024-02-07               0.47573    
as.factor(date)2024-02-21               0.44071    
as.factor(date)2024-03-06               0.77905    
as.factor(date)2024-03-22               0.26472    
as.factor(date)2024-04-04               0.56968    
as.factor(date)2024-04-16               0.30265    
treatmentMC1                            0.00575 ** 
treatmentMC2                            0.01642 *  
treatmentPC0                            0.27195    
treatmentPC1                            0.01314 *  
treatmentPC2                            0.09662 .  
as.factor(date)2024-01-12:treatmentMC1  0.24240    
as.factor(date)2024-01-23:treatmentMC1  0.04420 *  
as.factor(date)2024-02-07:treatmentMC1  0.16886    
as.factor(date)2024-02-21:treatmentMC1  0.18321    
as.factor(date)2024-03-06:treatmentMC1  0.13220    
as.factor(date)2024-03-22:treatmentMC1  0.07279 .  
as.factor(date)2024-04-04:treatmentMC1  0.16395    
as.factor(date)2024-04-16:treatmentMC1  0.19419    
as.factor(date)2024-01-12:treatmentMC2  0.97348    
as.factor(date)2024-01-23:treatmentMC2  0.68860    
as.factor(date)2024-02-07:treatmentMC2  0.45713    
as.factor(date)2024-02-21:treatmentMC2  0.72091    
as.factor(date)2024-03-06:treatmentMC2  0.84414    
as.factor(date)2024-03-22:treatmentMC2  0.83935    
as.factor(date)2024-04-04:treatmentMC2  0.43717    
as.factor(date)2024-04-16:treatmentMC2  0.87009    
as.factor(date)2024-01-12:treatmentPC0  0.45009    
as.factor(date)2024-01-23:treatmentPC0  0.59786    
as.factor(date)2024-02-07:treatmentPC0  0.34070    
as.factor(date)2024-02-21:treatmentPC0  0.64887    
as.factor(date)2024-03-06:treatmentPC0  0.48372    
as.factor(date)2024-03-22:treatmentPC0  0.61592    
as.factor(date)2024-04-04:treatmentPC0  0.70008    
as.factor(date)2024-04-16:treatmentPC0  0.25992    
as.factor(date)2024-01-12:treatmentPC1  0.70851    
as.factor(date)2024-01-23:treatmentPC1  0.25786    
as.factor(date)2024-02-07:treatmentPC1  0.22873    
as.factor(date)2024-02-21:treatmentPC1  0.23809    
as.factor(date)2024-03-06:treatmentPC1  0.36193    
as.factor(date)2024-03-22:treatmentPC1  0.06851 .  
as.factor(date)2024-04-04:treatmentPC1  0.26965    
as.factor(date)2024-04-16:treatmentPC1  0.46524    
as.factor(date)2024-01-12:treatmentPC2  0.85712    
as.factor(date)2024-01-23:treatmentPC2  0.33686    
as.factor(date)2024-02-07:treatmentPC2  0.98574    
as.factor(date)2024-02-21:treatmentPC2  0.62933    
as.factor(date)2024-03-06:treatmentPC2  0.74604    
as.factor(date)2024-03-22:treatmentPC2  0.50838    
as.factor(date)2024-04-04:treatmentPC2  0.59515    
as.factor(date)2024-04-16:treatmentPC2  0.99746    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:
 Formula: ~1 | site_block_id
        (Intercept)
StdDev:      0.2151

 Formula: ~1 | trt_site_block_id %in% site_block_id
        (Intercept) Residual
StdDev:      0.1246   0.2419


Number of Observations: 398
Number of Groups: 
                       site_block_id trt_site_block_id %in% site_block_id 
                                   8                                   48 

logLik     df    AIC    BIC 
-86.80     57 287.60 506.52 

plot(fitted(fit_lme_ph), residuals(fit_lme_ph),
     xlab = "Fitted Values",
     ylab = "Residuals",
     main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)

emmip(fit_lme_ph, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = df_model,
             aes(x = date, y = p_h, color = treatment),
             alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  theme_bw()

Warning: Removed 34 rows containing missing values or values outside the scale range
(`geom_point()`).

Left-Censored Mixed Model Testing (brms)

Data

response_vars <- c("nh4", "no3", "mg", "p")

lc_model_data <- df_combined %>%
  mutate(
    nh4_censored = ifelse(nh4 == "LOD", -1, 0),
    no3_censored = ifelse(no3 == "LOD", -1, 0),
    mg_censored = ifelse(mg == "LOD", -1, 0),
    p_censored = ifelse(p == "LOD", -1, 0),
    
    nh4 = ifelse(nh4 == "LOD", nh4_lod, as.numeric(nh4)),
    no3 = ifelse(no3 == "LOD", n03_lod, as.numeric(no3)),
    mg = ifelse(mg == "LOD", mg_lod, as.numeric(mg)),
    p = ifelse(p == "LOD", p_lod, as.numeric(p)),
    
    date = factor(date),
    treatment = factor(treatment),
    site = factor(site),
    block = factor(block),
    
    # Create combined grouping factors
    site_block_id = interaction(site, block, drop = TRUE),
    trt_site_block_id = interaction(treatment, site, block, drop = TRUE)
  )

Warning: There were 4 warnings in `mutate()`.
The first warning was:
ℹ In argument: `nh4 = ifelse(nh4 == "LOD", nh4_lod, as.numeric(nh4))`.
Caused by warning in `ifelse()`:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

Run and Save Models

fit_no3 <- brm(
  formula = bf(no3 | cens(no3_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
  data = lc_model_data,
  family = Gamma(link = "log"),
  chains = 4,
  warmup = 1000,
  iter = 5000,
  cores = parallel::detectCores(),
  seed = 12345,
  file = here::here("models", "fit_no3")
)

fit_nh4 <- brm(
  formula = bf(nh4 | cens(nh4_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
  data = lc_model_data,
  family = Gamma(link = "log"),
  chains = 4,
  warmup = 1000,
  iter = 5000,
  cores = parallel::detectCores(),
  seed = 12345,
  file = here::here("models", "fit_nh4")
)

fit_mg <- brm(
  formula = bf(mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
  data = lc_model_data,
  family = Gamma(link = "log"),
  chains = 4,
  warmup = 1000,
  iter = 5000,
  cores = parallel::detectCores(),
  seed = 12345,
  file = here::here("models", "fit_mg")
)

fit_p <- brm(
  formula = bf(p | cens(p_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
  data = lc_model_data,
  family = Gamma(link = "log"),
  chains = 4,
  warmup = 1000,
  iter = 5000,
  cores = parallel::detectCores(),
  seed = 12345,
  file = here::here("models", "fit_p")
)

fit_no3 <- read_rds(here::here("models", "fit_no3.rds"))
fit_nh4 <- read_rds(here::here("models", "fit_nh4.rds"))
fit_mg <- read_rds(here::here("models", "fit_mg.rds"))
fit_p <- read_rds(here::here("models", "fit_p.rds"))

NO3

summary(fit_no3)

Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.

Warning: There were 15861 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: no3 | cens(no3_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id) 
   Data: lc_model_data (Number of observations: 656) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~site_block_id (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.33      0.62     0.53     3.02 1.03       94      128

~trt_site_block_id (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.41      0.24     1.02     1.98 1.05       67      248

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                       6.72      1.26     4.43     9.32 1.62        7
date2024M01M12                 -1.56      1.16    -4.26     0.23 1.41        9
date2024M01M23                 -3.72      1.45    -6.79    -1.18 1.38        9
date2024M02M07                 -6.60      1.41    -9.62    -4.18 1.37        9
date2024M02M21               -194.86    127.24  -502.21   -20.14 1.14       35
date2024M03M06               -178.57    186.66  -679.90   -14.12 2.19        5
date2024M03M22                 -5.25      1.49    -8.13    -2.35 1.32       10
date2024M04M04               -254.08    182.26  -644.58   -22.22 1.11       27
date2024M04M16               -167.18    145.07  -580.69   -17.34 1.41        9
date2025M02M08                 -4.86      1.30    -7.65    -2.70 1.44        8
date2025M02M15                 -5.53      1.26    -8.07    -3.27 1.40        9
date2025M03M04                 -6.24      1.39    -9.18    -3.53 1.37        9
date2025M03M16                 -7.13      1.37   -10.06    -4.82 1.34       10
date2025M03M27               -346.88    245.23  -705.98   -21.22 1.71        6
date2025M04M11               -236.76    114.45  -422.86   -19.60 1.56        7
treatmentMC1                   -1.92      1.50    -5.04     0.68 1.37        9
treatmentMC2                   -0.79      1.63    -4.28     2.28 1.44        8
treatmentPC0                   -1.74      1.41    -4.51     0.88 1.42        9
treatmentPC1                   -0.84      1.55    -4.37     1.84 1.37        9
treatmentPC2                   -1.17      1.69    -4.29     2.49 1.43        9
date2024M01M12:treatmentMC1     1.50      1.47    -1.27     4.48 1.19       17
date2024M01M23:treatmentMC1     2.40      1.78    -1.07     5.88 1.22       14
date2024M02M07:treatmentMC1     2.92      1.80    -0.53     6.58 1.20       15
date2024M02M21:treatmentMC1  -117.50    222.02  -552.42   262.24 1.06       67
date2024M03M06:treatmentMC1  -143.62    273.77  -612.54   501.84 1.25       12
date2024M03M22:treatmentMC1    -0.51      1.99    -4.46     3.35 1.17       18
date2024M04M04:treatmentMC1   -94.38    287.42  -585.97   503.19 1.06       77
date2024M04M16:treatmentMC1  -141.80    251.98  -582.65   430.99 1.22       14
date2025M02M08:treatmentMC1     1.47      1.80    -1.85     5.39 1.20       15
date2025M02M15:treatmentMC1     2.54      1.70    -0.77     5.91 1.21       14
date2025M03M04:treatmentMC1     1.98      1.76    -1.43     5.64 1.18       18
date2025M03M16:treatmentMC1     1.34      1.95    -2.25     5.27 1.15       21
date2025M03M27:treatmentMC1    -3.68    311.52  -599.71   583.58 1.36        9
date2025M04M11:treatmentMC1  -122.96    244.67  -558.55   312.67 1.12       23
date2024M01M12:treatmentMC2     0.07      1.57    -2.78     3.47 1.24       13
date2024M01M23:treatmentMC2    -1.75      1.88    -5.18     2.01 1.25       12
date2024M02M07:treatmentMC2    -0.03      1.96    -3.70     3.98 1.22       14
date2024M02M21:treatmentMC2   187.66    127.30    12.93   493.52 1.14       35
date2024M03M06:treatmentMC2   170.50    186.79     5.28   672.38 2.21        5
date2024M03M22:treatmentMC2  -308.11    206.33  -689.71   -16.19 1.02      201
date2024M04M04:treatmentMC2   -76.98    270.95  -558.10   487.48 1.11       25
date2024M04M16:treatmentMC2   161.78    145.04    11.98   575.04 1.42        9
date2025M02M08:treatmentMC2     1.05      1.69    -2.26     4.48 1.28       11
date2025M02M15:treatmentMC2     0.80      1.65    -2.26     4.22 1.25       13
date2025M03M04:treatmentMC2     0.71      1.85    -2.70     4.68 1.23       14
date2025M03M16:treatmentMC2     0.84      1.88    -2.68     4.68 1.21       15
date2025M03M27:treatmentMC2   339.71    245.16    13.96   698.47 1.71        6
date2025M04M11:treatmentMC2  -125.94    223.27  -536.48   286.52 1.11       30
date2024M01M12:treatmentPC0     0.60      1.47    -1.99     3.59 1.22       13
date2024M01M23:treatmentPC0     1.69      1.81    -1.81     5.26 1.28       11
date2024M02M07:treatmentPC0  -306.17    199.86  -678.13   -12.13 1.02      277
date2024M02M21:treatmentPC0   189.64    127.30    14.95   496.61 1.14       35
date2024M03M06:treatmentPC0  -142.05    258.85  -630.23   388.00 1.26       11
date2024M03M22:treatmentPC0  -327.86    207.19  -686.88   -13.66 1.03      100
date2024M04M04:treatmentPC0   -93.27    265.58  -556.34   454.17 1.06       59
date2024M04M16:treatmentPC0  -174.47    262.01  -613.31   395.92 1.24       13
date2025M02M08:treatmentPC0     3.10      1.62     0.34     6.50 1.26       12
date2025M02M15:treatmentPC0     3.71      1.56     0.69     6.83 1.23       14
date2025M03M04:treatmentPC0     3.88      1.81     0.44     7.66 1.22       15
date2025M03M16:treatmentPC0     2.81      1.70    -0.39     6.35 1.22       15
date2025M03M27:treatmentPC0    -0.10    331.60  -587.85   590.65 1.51        7
date2025M04M11:treatmentPC0   232.20    114.61    14.41   417.56 1.57        7
date2024M01M12:treatmentPC1     0.20      1.39    -2.22     3.09 1.28       11
date2024M01M23:treatmentPC1    -1.63      1.85    -5.10     2.03 1.27       11
date2024M02M07:treatmentPC1     0.31      1.82    -2.92     4.01 1.18       16
date2024M02M21:treatmentPC1   187.80    127.19    13.05   494.90 1.14       35
date2024M03M06:treatmentPC1   172.35    186.73     7.62   673.57 2.18        5
date2024M03M22:treatmentPC1    -1.94      2.01    -5.90     2.03 1.15       21
date2024M04M04:treatmentPC1   -61.78    260.22  -531.84   469.68 1.06       57
date2024M04M16:treatmentPC1   160.68    145.03    10.59   573.66 1.41        9
date2025M02M08:treatmentPC1     1.17      1.73    -1.85     4.81 1.23       13
date2025M02M15:treatmentPC1     0.76      1.59    -2.14     3.92 1.22       14
date2025M03M04:treatmentPC1     0.77      1.86    -2.64     4.47 1.18       16
date2025M03M16:treatmentPC1     0.01      1.85    -3.51     3.78 1.14       21
date2025M03M27:treatmentPC1   342.26    245.26    17.22   701.38 1.71        6
date2025M04M11:treatmentPC1   233.82    114.59    16.59   418.52 1.56        7
date2024M01M12:treatmentPC2     0.84      1.78    -3.13     4.49 1.25       13
date2024M01M23:treatmentPC2     0.97      1.95    -3.35     4.88 1.24       14
date2024M02M07:treatmentPC2  -306.38    190.07  -669.51   -14.12 1.03      150
date2024M02M21:treatmentPC2  -155.88    238.12  -584.40   338.19 1.06       65
date2024M03M06:treatmentPC2   173.86    186.78    10.09   675.60 2.17        5
date2024M03M22:treatmentPC2     0.69      2.07    -3.90     4.62 1.24       15
date2024M04M04:treatmentPC2   250.07    182.36    17.75   640.94 1.11       27
date2024M04M16:treatmentPC2  -166.46    245.41  -593.25   357.19 1.19       17
date2025M02M08:treatmentPC2     1.33      1.83    -2.53     4.87 1.35       10
date2025M02M15:treatmentPC2     2.20      1.82    -1.61     5.66 1.38       10
date2025M03M04:treatmentPC2     1.07      2.01    -4.02     4.75 1.29       11
date2025M03M16:treatmentPC2     2.14      2.00    -2.38     6.12 1.30       12
date2025M03M27:treatmentPC2   342.51    245.55    17.33   701.86 1.71        6
date2025M04M11:treatmentPC2   231.90    114.55    14.88   419.06 1.56        7
                            Tail_ESS
Intercept                         22
date2024M01M12                    30
date2024M01M23                    22
date2024M02M07                    28
date2024M02M21                    23
date2024M03M06                    15
date2024M03M22                   105
date2024M04M04                    87
date2024M04M16                    16
date2025M02M08                    20
date2025M02M15                    45
date2025M03M04                    29
date2025M03M16                    41
date2025M03M27                    55
date2025M04M11                    22
treatmentMC1                      38
treatmentMC2                      25
treatmentPC0                      21
treatmentPC1                      28
treatmentPC2                      19
date2024M01M12:treatmentMC1       61
date2024M01M23:treatmentMC1       35
date2024M02M07:treatmentMC1       62
date2024M02M21:treatmentMC1      144
date2024M03M06:treatmentMC1       22
date2024M03M22:treatmentMC1      179
date2024M04M04:treatmentMC1      103
date2024M04M16:treatmentMC1       14
date2025M02M08:treatmentMC1       46
date2025M02M15:treatmentMC1      122
date2025M03M04:treatmentMC1       69
date2025M03M16:treatmentMC1       86
date2025M03M27:treatmentMC1       99
date2025M04M11:treatmentMC1      159
date2024M01M12:treatmentMC2       43
date2024M01M23:treatmentMC2       36
date2024M02M07:treatmentMC2       32
date2024M02M21:treatmentMC2       23
date2024M03M06:treatmentMC2       15
date2024M03M22:treatmentMC2      450
date2024M04M04:treatmentMC2      148
date2024M04M16:treatmentMC2       16
date2025M02M08:treatmentMC2       26
date2025M02M15:treatmentMC2       39
date2025M03M04:treatmentMC2       71
date2025M03M16:treatmentMC2       34
date2025M03M27:treatmentMC2       55
date2025M04M11:treatmentMC2       53
date2024M01M12:treatmentPC0       43
date2024M01M23:treatmentPC0       38
date2024M02M07:treatmentPC0      836
date2024M02M21:treatmentPC0       23
date2024M03M06:treatmentPC0       25
date2024M03M22:treatmentPC0      282
date2024M04M04:treatmentPC0      152
date2024M04M16:treatmentPC0       13
date2025M02M08:treatmentPC0       23
date2025M02M15:treatmentPC0       33
date2025M03M04:treatmentPC0       41
date2025M03M16:treatmentPC0       29
date2025M03M27:treatmentPC0       63
date2025M04M11:treatmentPC0       22
date2024M01M12:treatmentPC1       61
date2024M01M23:treatmentPC1       42
date2024M02M07:treatmentPC1      206
date2024M02M21:treatmentPC1       22
date2024M03M06:treatmentPC1       15
date2024M03M22:treatmentPC1       69
date2024M04M04:treatmentPC1      182
date2024M04M16:treatmentPC1       16
date2025M02M08:treatmentPC1       72
date2025M02M15:treatmentPC1       97
date2025M03M04:treatmentPC1       89
date2025M03M16:treatmentPC1      157
date2025M03M27:treatmentPC1       55
date2025M04M11:treatmentPC1       22
date2024M01M12:treatmentPC2       41
date2024M01M23:treatmentPC2       28
date2024M02M07:treatmentPC2      458
date2024M02M21:treatmentPC2       54
date2024M03M06:treatmentPC2       15
date2024M03M22:treatmentPC2       82
date2024M04M04:treatmentPC2       87
date2024M04M16:treatmentPC2       17
date2025M02M08:treatmentPC2       30
date2025M02M15:treatmentPC2       35
date2025M03M04:treatmentPC2       38
date2025M03M16:treatmentPC2       41
date2025M03M27:treatmentPC2       55
date2025M04M11:treatmentPC2       21

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.37      0.04     0.30     0.46 1.01      269      573

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

# neff_ratio(fit_no3)
# prior_summary(fit_no3)
# plot(fit, type = "trace", nvariables = 1)

pp_check(fit_no3) + xlim(0, 100)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 351 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_density()`).

pp_check(fit_no3, type = "dens_overlay", ndraws = 100) + xlim(0, 100)

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 3810 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_density()`).

# loo(fit_no3)

conditional_effects(fit_no3, effects = "treatment")

conditional_effects(fit_no3, effects = "date:treatment")

emm_treatments <- emmeans(fit_no3, ~ treatment, type = "response")

NOTE: Results may be misleading due to involvement in interactions

emm_treatments

 treatment response lower.HPD upper.HPD
 MC0       0.00e+00  2.22e-16      0.00
 MC1       0.00e+00  2.22e-16      0.00
 MC2       0.00e+00  2.22e-16      0.00
 PC0       0.00e+00  2.22e-16      0.00
 PC1       1.41e-08  2.22e-16      0.49
 PC2       0.00e+00  2.22e-16      0.00

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

pairs(emm_treatments)

 contrast     ratio lower.HPD upper.HPD
 MC0 / MC1 2.41e+18  2.22e-16  2.13e+50
 MC0 / MC2 0.00e+00  2.22e-16  7.13e+18
 MC0 / PC0 1.22e+19  2.22e-16  1.61e+50
 MC0 / PC1 0.00e+00  2.22e-16  0.00e+00
 MC0 / PC2 0.00e+00  2.22e-16  1.32e+15
 MC1 / MC2 0.00e+00  2.22e-16  0.00e+00
 MC1 / PC0 0.00e+00  2.22e-16  3.86e+33
 MC1 / PC1 0.00e+00  2.22e-16  0.00e+00
 MC1 / PC2 0.00e+00  2.22e-16  0.00e+00
 MC2 / PC0 1.10e+27  2.22e-16  2.29e+59
 MC2 / PC1 0.00e+00  2.22e-16  1.00e+00
 MC2 / PC2 0.00e+00  2.22e-16  4.56e+22
 PC0 / PC1 0.00e+00  2.22e-16  0.00e+00
 PC0 / PC2 0.00e+00  2.22e-16  1.00e+00
 PC1 / PC2 8.26e+18  4.07e-16  2.71e+38

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

emmip(fit_no3, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data,
             aes(x = date, y = no3, color = as.factor(no3_censored)),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  # ylim(0, 1000) +
  theme_bw()

Warning: Removed 64 rows containing missing values or values outside the scale range
(`geom_point()`).

NH4

summary(fit_nh4)

Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.

Warning: There were 15995 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: nh4 | cens(nh4_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id) 
   Data: lc_model_data (Number of observations: 663) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~site_block_id (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.51      0.41     0.01     1.56 1.05       47       88

~trt_site_block_id (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.86      0.36     1.20     2.74 1.34        9       21

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                       0.27      0.98    -1.45     2.35 1.23       15
date2024M01M12                 -1.64      1.15    -4.14     0.37 1.11       30
date2024M01M23               -110.65    142.37  -410.51    -6.52 1.68        6
date2024M02M07               -246.28    182.22  -671.07   -15.53 1.49        8
date2024M02M21               -449.79    193.63  -696.20   -48.32 1.45        8
date2024M03M06                 -0.92      1.32    -3.86     1.42 1.09       38
date2024M03M22                 -3.46      1.66    -6.67    -0.15 1.07       54
date2024M04M04                 -0.96      1.00    -3.11     0.85 1.05       57
date2024M04M16                 -3.16      1.14    -5.42    -0.94 1.08       45
date2025M02M08                 -1.89      1.07    -3.89     0.31 1.22       16
date2025M02M15                 -3.30      1.49    -6.33    -0.65 1.22       14
date2025M03M04                 -1.91      1.49    -4.39     1.92 1.22       15
date2025M03M16                 -3.03      1.18    -5.51    -0.84 1.12       25
date2025M03M27                 -1.25      1.08    -3.28     1.02 1.13       34
date2025M04M11                -47.63     61.30  -183.79    -3.44 1.87        6
treatmentMC1                    0.97      1.55    -1.71     4.46 1.18       16
treatmentMC2                   -0.41      1.30    -3.17     1.87 1.15       23
treatmentPC0                    0.59      1.22    -1.97     2.96 1.09       45
treatmentPC1                    0.28      1.52    -2.88     3.12 1.24       16
treatmentPC2                   -1.09      1.68    -4.40     2.27 1.17       19
date2024M01M12:treatmentMC1    -1.64      1.80    -5.12     2.10 1.15       21
date2024M01M23:treatmentMC1   106.70    142.45     2.75   405.97 1.69        6
date2024M02M07:treatmentMC1   -82.83    252.90  -563.34   404.13 1.24       13
date2024M02M21:treatmentMC1   149.14    266.31  -424.06   598.51 1.33       10
date2024M03M06:treatmentMC1  -414.26    193.97  -699.10   -32.68 1.25       12
date2024M03M22:treatmentMC1  -229.75    196.40  -672.44    -9.15 1.12       27
date2024M04M04:treatmentMC1    -2.13      1.69    -5.69     0.96 1.07       56
date2024M04M16:treatmentMC1  -300.57    197.17  -681.31   -14.27 1.01      168
date2025M02M08:treatmentMC1    -1.41      1.61    -4.79     1.51 1.09       34
date2025M02M15:treatmentMC1  -327.04    199.96  -680.54   -14.20 1.04       88
date2025M03M04:treatmentMC1    -1.13      2.10    -5.90     2.40 1.32       10
date2025M03M16:treatmentMC1  -338.93    205.82  -685.30   -16.92 1.13       30
date2025M03M27:treatmentMC1  -335.21    203.80  -673.01   -15.45 1.12       25
date2025M04M11:treatmentMC1  -302.47    220.82  -671.75   116.09 1.14       21
date2024M01M12:treatmentMC2     0.49      1.71    -2.93     3.98 1.15       31
date2024M01M23:treatmentMC2   107.64    142.51     3.97   408.47 1.68        6
date2024M02M07:treatmentMC2   -70.20    300.91  -601.10   545.05 1.46        8
date2024M02M21:treatmentMC2   107.52    279.08  -453.36   596.20 1.26       11
date2024M03M06:treatmentMC2  -355.19    210.25  -687.27   -16.92 1.12       23
date2024M03M22:treatmentMC2  -361.38    206.46  -685.03   -19.29 1.02      118
date2024M04M04:treatmentMC2    -1.02      1.67    -4.43     2.16 1.09       45
date2024M04M16:treatmentMC2  -343.27    199.20  -680.85   -18.64 1.05       89
date2025M02M08:treatmentMC2     0.23      1.46    -2.56     3.37 1.17       17
date2025M02M15:treatmentMC2  -315.23    189.98  -676.07   -17.03 1.02      147
date2025M03M04:treatmentMC2    -0.09      2.06    -4.96     3.70 1.19       18
date2025M03M16:treatmentMC2  -330.48    198.16  -684.94   -13.04 1.01      139
date2025M03M27:treatmentMC2     0.37      1.50    -2.47     3.26 1.14       25
date2025M04M11:treatmentMC2    46.89     61.36     1.55   182.38 1.85        6
date2024M01M12:treatmentPC0  -340.87    203.12  -690.13   -21.12 1.10       31
date2024M01M23:treatmentPC0  -214.39    251.67  -641.08   320.56 1.26       12
date2024M02M07:treatmentPC0   -94.06    268.29  -572.14   456.48 1.22       13
date2024M02M21:treatmentPC0   115.00    263.95  -383.78   622.94 1.18       16
date2024M03M06:treatmentPC0  -347.83    196.30  -680.71   -20.69 1.02      195
date2024M03M22:treatmentPC0  -317.88    221.06  -696.89   -11.40 1.09       36
date2024M04M04:treatmentPC0    -3.11      1.69    -6.55     0.26 1.04       84
date2024M04M16:treatmentPC0    -0.64      1.79    -3.98     2.95 1.11       42
date2025M02M08:treatmentPC0     0.11      1.70    -2.74     3.59 1.20       16
date2025M02M15:treatmentPC0    -0.35      2.15    -4.26     4.13 1.26       13
date2025M03M04:treatmentPC0    -0.35      2.07    -5.13     3.50 1.23       13
date2025M03M16:treatmentPC0    -1.21      1.85    -4.48     2.81 1.08       47
date2025M03M27:treatmentPC0    -0.79      1.72    -4.07     2.67 1.11       33
date2025M04M11:treatmentPC0    46.03     61.11     1.16   181.12 1.85        6
date2024M01M12:treatmentPC1    -2.28      1.79    -5.52     1.72 1.07       50
date2024M01M23:treatmentPC1   108.06    142.48     4.36   408.26 1.69        6
date2024M02M07:treatmentPC1   -67.53    276.98  -558.40   575.11 1.32       11
date2024M02M21:treatmentPC1   445.63    193.65    44.10   692.41 1.45        8
date2024M03M06:treatmentPC1  -330.59    194.33  -674.70   -17.14 1.01       68
date2024M03M22:treatmentPC1  -300.04    195.33  -677.19   -17.34 1.06      103
date2024M04M04:treatmentPC1    -0.99      1.77    -4.44     2.52 1.05       64
date2024M04M16:treatmentPC1     1.52      1.85    -2.07     4.99 1.08       49
date2025M02M08:treatmentPC1    -0.49      1.86    -3.90     3.09 1.13       30
date2025M02M15:treatmentPC1     0.26      2.01    -3.25     4.40 1.20       16
date2025M03M04:treatmentPC1     0.93      2.04    -3.94     4.70 1.19       18
date2025M03M16:treatmentPC1     0.96      1.89    -2.70     4.86 1.10       30
date2025M03M27:treatmentPC1     0.90      1.80    -2.56     4.34 1.10       36
date2025M04M11:treatmentPC1    48.23     61.39     3.57   184.31 1.92        6
date2024M01M12:treatmentPC2    -1.58      1.85    -5.50     1.81 1.12       24
date2024M01M23:treatmentPC2   109.03    142.38     3.50   409.48 1.68        6
date2024M02M07:treatmentPC2   243.45    182.27    13.00   668.96 1.49        8
date2024M02M21:treatmentPC2    76.69    271.22  -564.20   547.48 1.23       13
date2024M03M06:treatmentPC2  -348.29    204.86  -685.51   -17.95 1.02      103
date2024M03M22:treatmentPC2  -278.55    190.88  -674.71   -14.59 1.10       38
date2024M04M04:treatmentPC2    -1.89      2.03    -5.83     2.15 1.08       54
date2024M04M16:treatmentPC2  -314.52    188.60  -668.68   -16.19 1.03      153
date2025M02M08:treatmentPC2    -0.28      1.92    -4.01     3.25 1.21       17
date2025M02M15:treatmentPC2  -341.83    200.21  -684.05   -13.71 1.03       97
date2025M03M04:treatmentPC2    -0.81      2.10    -5.94     2.73 1.13       26
date2025M03M16:treatmentPC2     2.10      1.84    -1.91     5.68 1.23       13
date2025M03M27:treatmentPC2     0.60      1.89    -3.14     4.29 1.18       19
date2025M04M11:treatmentPC2    48.59     61.07     3.87   184.03 1.87        6
                            Tail_ESS
Intercept                         79
date2024M01M12                   120
date2024M01M23                    21
date2024M02M07                    15
date2024M02M21                    22
date2024M03M06                   125
date2024M03M22                   100
date2024M04M04                    96
date2024M04M16                    88
date2025M02M08                    52
date2025M02M15                    67
date2025M03M04                    15
date2025M03M16                   118
date2025M03M27                   120
date2025M04M11                    20
treatmentMC1                      50
treatmentMC2                      44
treatmentPC0                      88
treatmentPC1                      71
treatmentPC2                      36
date2024M01M12:treatmentMC1      134
date2024M01M23:treatmentMC1       20
date2024M02M07:treatmentMC1       65
date2024M02M21:treatmentMC1       71
date2024M03M06:treatmentMC1       98
date2024M03M22:treatmentMC1      164
date2024M04M04:treatmentMC1      105
date2024M04M16:treatmentMC1      387
date2025M02M08:treatmentMC1       97
date2025M02M15:treatmentMC1      358
date2025M03M04:treatmentMC1       18
date2025M03M16:treatmentMC1      229
date2025M03M27:treatmentMC1      321
date2025M04M11:treatmentMC1       36
date2024M01M12:treatmentMC2       63
date2024M01M23:treatmentMC2       20
date2024M02M07:treatmentMC2       14
date2024M02M21:treatmentMC2       47
date2024M03M06:treatmentMC2      296
date2024M03M22:treatmentMC2      361
date2024M04M04:treatmentMC2      118
date2024M04M16:treatmentMC2      461
date2025M02M08:treatmentMC2       89
date2025M02M15:treatmentMC2      356
date2025M03M04:treatmentMC2       41
date2025M03M16:treatmentMC2      342
date2025M03M27:treatmentMC2      124
date2025M04M11:treatmentMC2       19
date2024M01M12:treatmentPC0      178
date2024M01M23:treatmentPC0       31
date2024M02M07:treatmentPC0       24
date2024M02M21:treatmentPC0       58
date2024M03M06:treatmentPC0      480
date2024M03M22:treatmentPC0      227
date2024M04M04:treatmentPC0      137
date2024M04M16:treatmentPC0       96
date2025M02M08:treatmentPC0      108
date2025M02M15:treatmentPC0       34
date2025M03M04:treatmentPC0       13
date2025M03M16:treatmentPC0      122
date2025M03M27:treatmentPC0      118
date2025M04M11:treatmentPC0       20
date2024M01M12:treatmentPC1      118
date2024M01M23:treatmentPC1       21
date2024M02M07:treatmentPC1       12
date2024M02M21:treatmentPC1       22
date2024M03M06:treatmentPC1      247
date2024M03M22:treatmentPC1      300
date2024M04M04:treatmentPC1      138
date2024M04M16:treatmentPC1       69
date2025M02M08:treatmentPC1      125
date2025M02M15:treatmentPC1       67
date2025M03M04:treatmentPC1       20
date2025M03M16:treatmentPC1      130
date2025M03M27:treatmentPC1       46
date2025M04M11:treatmentPC1       18
date2024M01M12:treatmentPC2      107
date2024M01M23:treatmentPC2       20
date2024M02M07:treatmentPC2       15
date2024M02M21:treatmentPC2       40
date2024M03M06:treatmentPC2      489
date2024M03M22:treatmentPC2      477
date2024M04M04:treatmentPC2      120
date2024M04M16:treatmentPC2      390
date2025M02M08:treatmentPC2       52
date2025M02M15:treatmentPC2      400
date2025M03M04:treatmentPC2       29
date2025M03M16:treatmentPC2       30
date2025M03M27:treatmentPC2       73
date2025M04M11:treatmentPC2       18

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.30      0.05     0.21     0.40 1.08       47      212

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

# neff_ratio(fit_nh4)
# prior_summary(fit_nh4)
# plot(fit, type = "trace", nvariables = 1)

pp_check(fit_nh4) + xlim(0, 10)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 62 rows containing non-finite outside the scale range
(`stat_density()`).

pp_check(fit_nh4, type = "dens_overlay", ndraws = 100) + xlim(0, 10)

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 726 rows containing non-finite outside the scale range
(`stat_density()`).

# loo(fit_nh4)

conditional_effects(fit_nh4, effects = "treatment")

conditional_effects(fit_nh4, effects = "date:treatment")

emm_treatments <- emmeans(fit_nh4, ~ treatment, type = "response")

NOTE: Results may be misleading due to involvement in interactions

emm_treatments

 treatment   response  lower.HPD  upper.HPD
 MC0       2.2204e-16 2.2204e-16 2.1125e-13
 MC1       2.2204e-16 2.2204e-16 2.2200e-16
 MC2       2.2204e-16 2.2204e-16 2.2200e-16
 PC0       2.2204e-16 2.2204e-16 2.2200e-16
 PC1       2.2204e-16 2.2204e-16 7.5358e-13
 PC2       2.2204e-16 2.2204e-16 2.2200e-16

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

pairs(emm_treatments)

 contrast     ratio lower.HPD upper.HPD
 MC0 / MC1 3.18e+60    148877  1.60e+86
 MC0 / MC2 1.10e+43         0  2.11e+76
 MC0 / PC0 1.74e+34         0  8.32e+56
 MC0 / PC1 6.70e+01         0  2.13e+25
 MC0 / PC2 5.49e+23         0  3.70e+48
 MC1 / MC2 0.00e+00         0  2.32e+22
 MC1 / PC0 0.00e+00         0  1.82e+04
 MC1 / PC1 0.00e+00         0  0.00e+00
 MC1 / PC2 0.00e+00         0  0.00e+00
 MC2 / PC0 0.00e+00         0  1.66e+30
 MC2 / PC1 0.00e+00         0  0.00e+00
 MC2 / PC2 0.00e+00         0  1.14e+14
 PC0 / PC1 0.00e+00         0  0.00e+00
 PC0 / PC2 0.00e+00         0  5.06e+24
 PC1 / PC2 5.85e+21         0  9.11e+47

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

plot(emm_treatments) + theme_bw()

emmip(fit_nh4, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data,
             aes(x = date, y = nh4, color = as.factor(mg_censored)),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  # ylim(0, 10) +
  theme_bw()

Warning: Removed 57 rows containing missing values or values outside the scale range
(`geom_point()`).

Mg

summary(fit_mg)

Warning: There were 149 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id) 
   Data: lc_model_data (Number of observations: 633) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~site_block_id (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.18     0.17     0.87 1.00     2032      743

~trt_site_block_id (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.38      0.05     0.29     0.50 1.00     6068     9749

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                       3.90      0.30     3.29     4.45 1.00      869
date2024M01M12                  0.20      0.24    -0.28     0.68 1.00     2141
date2024M01M23                  0.25      0.24    -0.23     0.71 1.00     1933
date2024M02M07                 -0.01      0.24    -0.48     0.46 1.00     2122
date2024M02M21                  0.12      0.25    -0.36     0.62 1.00     2153
date2024M03M06                  0.19      0.25    -0.30     0.66 1.00     1806
date2024M03M22                  0.38      0.24    -0.09     0.86 1.00     1848
date2024M04M04                  0.48      0.26    -0.02     0.98 1.00     2203
date2024M04M16                  0.34      0.24    -0.14     0.81 1.00     2128
date2025M02M08                 -0.66      0.24    -1.14    -0.19 1.00     1980
date2025M02M15                 -0.48      0.25    -0.97     0.00 1.00     2085
date2025M03M04                 -0.35      0.26    -0.85     0.16 1.00     2299
date2025M03M16                 -0.14      0.24    -0.61     0.33 1.00     2104
date2025M03M27                  0.03      0.26    -0.47     0.53 1.00     2299
date2025M04M11                  0.21      0.27    -0.31     0.75 1.00     1922
treatmentMC1                    0.86      0.31     0.25     1.47 1.00     2513
treatmentMC2                    0.58      0.32    -0.04     1.22 1.00     2455
treatmentPC0                    0.01      0.32    -0.62     0.62 1.00     2577
treatmentPC1                    0.48      0.31    -0.12     1.09 1.00     2187
treatmentPC2                    0.10      0.35    -0.57     0.79 1.00     1728
date2024M01M12:treatmentMC1    -0.27      0.33    -0.92     0.38 1.00     2988
date2024M01M23:treatmentMC1    -0.32      0.33    -0.96     0.33 1.00     2830
date2024M02M07:treatmentMC1    -0.55      0.32    -1.20     0.07 1.00     2961
date2024M02M21:treatmentMC1    -0.73      0.33    -1.39    -0.09 1.00     2731
date2024M03M06:treatmentMC1    -0.81      0.33    -1.47    -0.15 1.00     2874
date2024M03M22:treatmentMC1    -0.80      0.33    -1.44    -0.15 1.00     2825
date2024M04M04:treatmentMC1    -0.92      0.35    -1.60    -0.23 1.00     2874
date2024M04M16:treatmentMC1    -1.01      0.33    -1.66    -0.37 1.00     2227
date2025M02M08:treatmentMC1     0.01      0.34    -0.67     0.66 1.00     2500
date2025M02M15:treatmentMC1    -0.07      0.33    -0.72     0.57 1.00     2770
date2025M03M04:treatmentMC1    -0.10      0.34    -0.78     0.56 1.00     3144
date2025M03M16:treatmentMC1    -0.47      0.32    -1.12     0.16 1.00     2787
date2025M03M27:treatmentMC1    -0.50      0.34    -1.19     0.16 1.00     2750
date2025M04M11:treatmentMC1    -0.40      0.36    -1.11     0.29 1.00     3045
date2024M01M12:treatmentMC2    -0.04      0.33    -0.69     0.62 1.00     2875
date2024M01M23:treatmentMC2     0.06      0.33    -0.59     0.70 1.00     2775
date2024M02M07:treatmentMC2    -0.19      0.33    -0.84     0.45 1.00     2798
date2024M02M21:treatmentMC2    -0.14      0.34    -0.82     0.51 1.00     2710
date2024M03M06:treatmentMC2    -0.10      0.34    -0.76     0.57 1.00     2825
date2024M03M22:treatmentMC2    -0.11      0.33    -0.77     0.53 1.00     2511
date2024M04M04:treatmentMC2    -0.19      0.34    -0.86     0.48 1.00     2886
date2024M04M16:treatmentMC2    -0.37      0.34    -1.04     0.29 1.00     2516
date2025M02M08:treatmentMC2     0.19      0.35    -0.49     0.86 1.00     2463
date2025M02M15:treatmentMC2     0.66      0.34    -0.01     1.31 1.00     2887
date2025M03M04:treatmentMC2     0.55      0.35    -0.15     1.23 1.00     3094
date2025M03M16:treatmentMC2     0.21      0.33    -0.44     0.85 1.00     2700
date2025M03M27:treatmentMC2     0.12      0.34    -0.56     0.79 1.00     3062
date2025M04M11:treatmentMC2     0.12      0.39    -0.64     0.88 1.00     3304
date2024M01M12:treatmentPC0    -0.31      0.33    -0.95     0.35 1.00     2992
date2024M01M23:treatmentPC0    -0.56      0.32    -1.19     0.07 1.00     2745
date2024M02M07:treatmentPC0    -0.40      0.33    -1.05     0.23 1.00     2922
date2024M02M21:treatmentPC0    -0.37      0.33    -1.02     0.28 1.00     3116
date2024M03M06:treatmentPC0    -0.36      0.33    -1.00     0.30 1.00     2483
date2024M03M22:treatmentPC0    -0.54      0.33    -1.19     0.11 1.00     2706
date2024M04M04:treatmentPC0    -0.59      0.34    -1.25     0.09 1.00     2953
date2024M04M16:treatmentPC0    -0.51      0.33    -1.17     0.13 1.00     3257
date2025M02M08:treatmentPC0    -0.59      0.39    -1.35     0.17 1.00     3075
date2025M02M15:treatmentPC0    -0.55      0.33    -1.21     0.10 1.00     2969
date2025M03M04:treatmentPC0    -0.46      0.36    -1.17     0.26 1.00     3014
date2025M03M16:treatmentPC0    -0.81      0.34    -1.47    -0.14 1.00     3140
date2025M03M27:treatmentPC0    -0.67      0.38    -1.39     0.08 1.00     3335
date2025M04M11:treatmentPC0    -0.53      0.37    -1.25     0.19 1.00     2247
date2024M01M12:treatmentPC1    -0.46      0.33    -1.11     0.19 1.00     2995
date2024M01M23:treatmentPC1    -0.60      0.33    -1.23     0.05 1.00     2892
date2024M02M07:treatmentPC1    -0.91      0.32    -1.54    -0.29 1.00     2752
date2024M02M21:treatmentPC1    -0.84      0.33    -1.48    -0.21 1.00     2476
date2024M03M06:treatmentPC1    -0.73      0.33    -1.37    -0.08 1.00     2880
date2024M03M22:treatmentPC1    -0.82      0.33    -1.46    -0.18 1.00     2772
date2024M04M04:treatmentPC1    -0.95      0.34    -1.61    -0.28 1.00     2979
date2024M04M16:treatmentPC1    -1.03      0.34    -1.70    -0.37 1.00     2945
date2025M02M08:treatmentPC1    -1.01      0.33    -1.66    -0.36 1.00     2785
date2025M02M15:treatmentPC1    -1.05      0.34    -1.72    -0.39 1.00     2882
date2025M03M04:treatmentPC1    -1.12      0.34    -1.80    -0.45 1.00     3182
date2025M03M16:treatmentPC1    -1.05      0.34    -1.71    -0.39 1.00     3087
date2025M03M27:treatmentPC1    -1.04      0.35    -1.72    -0.34 1.00     3212
date2025M04M11:treatmentPC1    -1.23      0.40    -2.00    -0.44 1.00     3125
date2024M01M12:treatmentPC2    -0.29      0.36    -1.00     0.43 1.00     2778
date2024M01M23:treatmentPC2    -0.55      0.36    -1.25     0.14 1.00     2333
date2024M02M07:treatmentPC2    -0.78      0.36    -1.47    -0.07 1.00     2159
date2024M02M21:treatmentPC2    -0.61      0.36    -1.31     0.11 1.00     2491
date2024M03M06:treatmentPC2    -0.58      0.37    -1.30     0.13 1.00     1856
date2024M03M22:treatmentPC2    -0.55      0.36    -1.26     0.15 1.00     2027
date2024M04M04:treatmentPC2    -0.62      0.37    -1.34     0.12 1.00     2706
date2024M04M16:treatmentPC2    -0.68      0.36    -1.39     0.02 1.00     3141
date2025M02M08:treatmentPC2    -0.65      0.37    -1.37     0.07 1.00     2916
date2025M02M15:treatmentPC2    -0.36      0.36    -1.07     0.35 1.00     2638
date2025M03M04:treatmentPC2    -0.22      0.38    -0.96     0.52 1.00     2541
date2025M03M16:treatmentPC2    -0.40      0.36    -1.11     0.30 1.00     2830
date2025M03M27:treatmentPC2    -0.38      0.38    -1.14     0.36 1.00     2334
date2025M04M11:treatmentPC2    -0.56      0.43    -1.38     0.30 1.00     2570
                            Tail_ESS
Intercept                        489
date2024M01M12                  4611
date2024M01M23                  4161
date2024M02M07                  4299
date2024M02M21                  5113
date2024M03M06                  5267
date2024M03M22                  5165
date2024M04M04                  5257
date2024M04M16                  5252
date2025M02M08                  4562
date2025M02M15                  5271
date2025M03M04                  5281
date2025M03M16                  5087
date2025M03M27                  5755
date2025M04M11                  5094
treatmentMC1                    5223
treatmentMC2                    6142
treatmentPC0                    6150
treatmentPC1                    5902
treatmentPC2                    1876
date2024M01M12:treatmentMC1     7291
date2024M01M23:treatmentMC1     6716
date2024M02M07:treatmentMC1     6611
date2024M02M21:treatmentMC1     5988
date2024M03M06:treatmentMC1     6631
date2024M03M22:treatmentMC1     6946
date2024M04M04:treatmentMC1     6233
date2024M04M16:treatmentMC1     5203
date2025M02M08:treatmentMC1     6029
date2025M02M15:treatmentMC1     5923
date2025M03M04:treatmentMC1     6228
date2025M03M16:treatmentMC1     5853
date2025M03M27:treatmentMC1     7492
date2025M04M11:treatmentMC1     7125
date2024M01M12:treatmentMC2     7081
date2024M01M23:treatmentMC2     6424
date2024M02M07:treatmentMC2     6043
date2024M02M21:treatmentMC2     6387
date2024M03M06:treatmentMC2     7069
date2024M03M22:treatmentMC2     5811
date2024M04M04:treatmentMC2     6399
date2024M04M16:treatmentMC2     7370
date2025M02M08:treatmentMC2     6477
date2025M02M15:treatmentMC2     6186
date2025M03M04:treatmentMC2     6360
date2025M03M16:treatmentMC2     6414
date2025M03M27:treatmentMC2     7132
date2025M04M11:treatmentMC2     6804
date2024M01M12:treatmentPC0     6866
date2024M01M23:treatmentPC0     7129
date2024M02M07:treatmentPC0     6789
date2024M02M21:treatmentPC0     6619
date2024M03M06:treatmentPC0     5943
date2024M03M22:treatmentPC0     7353
date2024M04M04:treatmentPC0     6545
date2024M04M16:treatmentPC0     7638
date2025M02M08:treatmentPC0     8050
date2025M02M15:treatmentPC0     6781
date2025M03M04:treatmentPC0     7508
date2025M03M16:treatmentPC0     6850
date2025M03M27:treatmentPC0     7941
date2025M04M11:treatmentPC0     5367
date2024M01M12:treatmentPC1     6645
date2024M01M23:treatmentPC1     6957
date2024M02M07:treatmentPC1     7158
date2024M02M21:treatmentPC1     6270
date2024M03M06:treatmentPC1     7322
date2024M03M22:treatmentPC1     6564
date2024M04M04:treatmentPC1     6755
date2024M04M16:treatmentPC1     7499
date2025M02M08:treatmentPC1     6374
date2025M02M15:treatmentPC1     7465
date2025M03M04:treatmentPC1     7572
date2025M03M16:treatmentPC1     7202
date2025M03M27:treatmentPC1     7874
date2025M04M11:treatmentPC1     8764
date2024M01M12:treatmentPC2     6428
date2024M01M23:treatmentPC2     6396
date2024M02M07:treatmentPC2     6143
date2024M02M21:treatmentPC2     6478
date2024M03M06:treatmentPC2     3113
date2024M03M22:treatmentPC2     4921
date2024M04M04:treatmentPC2     7305
date2024M04M16:treatmentPC2     6883
date2025M02M08:treatmentPC2     7278
date2025M02M15:treatmentPC2     7244
date2025M03M04:treatmentPC2     6715
date2025M03M16:treatmentPC2     6105
date2025M03M27:treatmentPC2     3233
date2025M04M11:treatmentPC2     7080

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     5.78      0.35     5.11     6.48 1.00    13355    11217

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

# neff_ratio(fit_mg)
# prior_summary(fit_mg)
# plot(fit, type = "trace")

pp_check(fit_mg)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

Warning: Censored responses are not shown in 'pp_check'.

pp_check(fit_mg, type = "dens_overlay", ndraws = 100)

Warning: Censored responses are not shown in 'pp_check'.

# loo(fit_mg)

conditional_effects(fit_mg, effects = "treatment")

conditional_effects(fit_mg, effects = "date:treatment")

emm_treatments <- emmeans(fit_mg, ~ treatment, type = "response")

NOTE: Results may be misleading due to involvement in interactions

emm_treatments

 treatment response lower.HPD upper.HPD
 MC0           51.7      31.1      77.2
 MC1           77.3      46.0     114.8
 MC2           97.6      59.8     145.5
 PC0           32.2      19.8      48.3
 PC1           35.7      21.5      52.9
 PC2           35.4      21.2      52.9

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

pairs(emm_treatments)

 contrast  ratio lower.HPD upper.HPD
 MC0 / MC1 0.670     0.433     0.967
 MC0 / MC2 0.530     0.336     0.755
 MC0 / PC0 1.608     1.039     2.333
 MC0 / PC1 1.447     0.930     2.082
 MC0 / PC2 1.462     0.932     2.091
 MC1 / MC2 0.791     0.512     1.130
 MC1 / PC0 2.407     1.534     3.463
 MC1 / PC1 2.161     1.391     3.082
 MC1 / PC2 2.176     1.375     3.086
 MC2 / PC0 3.036     1.931     4.361
 MC2 / PC1 2.728     1.766     3.907
 MC2 / PC2 2.757     1.785     3.921
 PC0 / PC1 0.901     0.581     1.297
 PC0 / PC2 0.906     0.578     1.294
 PC1 / PC2 1.006     0.652     1.443

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

plot(emm_treatments) + theme_bw()

emmip(fit_mg, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data,
             aes(x = date, y = mg, color = as.factor(mg_censored)),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  theme_bw()

Warning: Removed 87 rows containing missing values or values outside the scale range
(`geom_point()`).

fit_glmer_gamma_mg <- glmer(
  mg ~ date * treatment +
    (1 | site_block_id) +
    (1 | trt_site_block_id),
  data = lc_model_data %>% filter(study_year == 1),
  family = Gamma(link = "log"),
  # control = glmerControl(optimizer = "bobyqa"), # default optimizer did not converge
  na.action = na.exclude
)

Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
failure to converge in 10000 evaluations

Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0499573 (tol = 0.002, component 1)

summary(fit_glmer_gamma_mg)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
   Data: lc_model_data %>% filter(study_year == 1)

     AIC      BIC   logLik deviance df.resid 
  3664.3   3891.8  -1775.1   3550.3      343 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8607 -0.5774 -0.0216  0.4681  3.8763 

Random effects:
 Groups            Name        Variance Std.Dev.
 trt_site_block_id (Intercept) 0.04994  0.2235  
 site_block_id     (Intercept) 0.02173  0.1474  
 Residual                      0.12052  0.3472  
Number of obs: 400, groups:  trt_site_block_id, 48; site_block_id, 8

Fixed effects:
                            Estimate Std. Error t value Pr(>|z|)    
(Intercept)                  3.91105    0.25136  15.560  < 2e-16 ***
date2024-01-12               0.21160    0.20181   1.049 0.294392    
date2024-01-23               0.22794    0.19553   1.166 0.243703    
date2024-02-07              -0.03191    0.19650  -0.162 0.870998    
date2024-02-21               0.07750    0.20333   0.381 0.703102    
date2024-03-06               0.19660    0.20321   0.967 0.333314    
date2024-03-22               0.39571    0.20300   1.949 0.051257 .  
date2024-04-04               0.46917    0.21143   2.219 0.026480 *  
date2024-04-16               0.28751    0.20094   1.431 0.152480    
treatmentMC1                 0.81928    0.29224   2.803 0.005056 ** 
treatmentMC2                 0.54335    0.29847   1.820 0.068692 .  
treatmentPC0                -0.08276    0.29335  -0.282 0.777855    
treatmentPC1                 0.37059    0.29228   1.268 0.204831    
treatmentPC2                -0.14076    0.31568  -0.446 0.655674    
date2024-01-12:treatmentMC1 -0.27861    0.27063  -1.029 0.303254    
date2024-01-23:treatmentMC1 -0.27466    0.26672  -1.030 0.303118    
date2024-02-07:treatmentMC1 -0.54121    0.26338  -2.055 0.039889 *  
date2024-02-21:treatmentMC1 -0.69731    0.26914  -2.591 0.009572 ** 
date2024-03-06:treatmentMC1 -0.81882    0.26931  -3.040 0.002362 ** 
date2024-03-22:treatmentMC1 -0.82064    0.26925  -3.048 0.002305 ** 
date2024-04-04:treatmentMC1 -0.87750    0.28620  -3.066 0.002169 ** 
date2024-04-16:treatmentMC1 -0.97674    0.26785  -3.647 0.000266 ***
date2024-01-12:treatmentMC2 -0.11081    0.27352  -0.405 0.685374    
date2024-01-23:treatmentMC2  0.09168    0.26912   0.341 0.733342    
date2024-02-07:treatmentMC2 -0.17715    0.27082  -0.654 0.513026    
date2024-02-21:treatmentMC2 -0.10351    0.27572  -0.375 0.707338    
date2024-03-06:treatmentMC2 -0.10454    0.27602  -0.379 0.704894    
date2024-03-22:treatmentMC2 -0.12819    0.27576  -0.465 0.642026    
date2024-04-04:treatmentMC2 -0.17171    0.28198  -0.609 0.542548    
date2024-04-16:treatmentMC2 -0.29992    0.27906  -1.075 0.282495    
date2024-01-12:treatmentPC0 -0.26652    0.26820  -0.994 0.320350    
date2024-01-23:treatmentPC0 -0.49111    0.26410  -1.860 0.062950 .  
date2024-02-07:treatmentPC0 -0.34041    0.26482  -1.285 0.198634    
date2024-02-21:treatmentPC0 -0.28518    0.26958  -1.058 0.290120    
date2024-03-06:treatmentPC0 -0.32802    0.26941  -1.218 0.223397    
date2024-03-22:treatmentPC0 -0.50855    0.26957  -1.887 0.059222 .  
date2024-04-04:treatmentPC0 -0.51908    0.27641  -1.878 0.060388 .  
date2024-04-16:treatmentPC0 -0.40947    0.26879  -1.523 0.127658    
date2024-01-12:treatmentPC1 -0.44969    0.27120  -1.658 0.097280 .  
date2024-01-23:treatmentPC1 -0.54340    0.26740  -2.032 0.042137 *  
date2024-02-07:treatmentPC1 -0.83036    0.26432  -3.142 0.001681 ** 
date2024-02-21:treatmentPC1 -0.71730    0.26926  -2.664 0.007724 ** 
date2024-03-06:treatmentPC1 -0.65606    0.26952  -2.434 0.014924 *  
date2024-03-22:treatmentPC1 -0.73684    0.26996  -2.729 0.006345 ** 
date2024-04-04:treatmentPC1 -0.81262    0.28175  -2.884 0.003925 ** 
date2024-04-16:treatmentPC1 -0.91260    0.27643  -3.301 0.000962 ***
date2024-01-12:treatmentPC2 -0.08284    0.29625  -0.280 0.779770    
date2024-01-23:treatmentPC2 -0.34560    0.28816  -1.199 0.230409    
date2024-02-07:treatmentPC2 -0.52254    0.29101  -1.796 0.072560 .  
date2024-02-21:treatmentPC2 -0.35759    0.29502  -1.212 0.225484    
date2024-03-06:treatmentPC2 -0.37510    0.29881  -1.255 0.209362    
date2024-03-22:treatmentPC2 -0.36350    0.29492  -1.233 0.217737    
date2024-04-04:treatmentPC2 -0.38965    0.30512  -1.277 0.201590    
date2024-04-16:treatmentPC2 -0.44361    0.29343  -1.512 0.130582    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 54 > 12.
Use print(value, correlation=TRUE)  or
    vcov(value)        if you need it

optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations)
Model failed to converge with max|grad| = 0.0499573 (tol = 0.002, component 1)
failure to converge in 10000 evaluations

S(fit_mg)

Warning: There were 149 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id) 
   Data: lc_model_data (Number of observations: 633) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~site_block_id (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.18     0.17     0.87 1.00     2032      743

~trt_site_block_id (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.38      0.05     0.29     0.50 1.00     6068     9749

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                       3.90      0.30     3.29     4.45 1.00      869
date2024M01M12                  0.20      0.24    -0.28     0.68 1.00     2141
date2024M01M23                  0.25      0.24    -0.23     0.71 1.00     1933
date2024M02M07                 -0.01      0.24    -0.48     0.46 1.00     2122
date2024M02M21                  0.12      0.25    -0.36     0.62 1.00     2153
date2024M03M06                  0.19      0.25    -0.30     0.66 1.00     1806
date2024M03M22                  0.38      0.24    -0.09     0.86 1.00     1848
date2024M04M04                  0.48      0.26    -0.02     0.98 1.00     2203
date2024M04M16                  0.34      0.24    -0.14     0.81 1.00     2128
date2025M02M08                 -0.66      0.24    -1.14    -0.19 1.00     1980
date2025M02M15                 -0.48      0.25    -0.97     0.00 1.00     2085
date2025M03M04                 -0.35      0.26    -0.85     0.16 1.00     2299
date2025M03M16                 -0.14      0.24    -0.61     0.33 1.00     2104
date2025M03M27                  0.03      0.26    -0.47     0.53 1.00     2299
date2025M04M11                  0.21      0.27    -0.31     0.75 1.00     1922
treatmentMC1                    0.86      0.31     0.25     1.47 1.00     2513
treatmentMC2                    0.58      0.32    -0.04     1.22 1.00     2455
treatmentPC0                    0.01      0.32    -0.62     0.62 1.00     2577
treatmentPC1                    0.48      0.31    -0.12     1.09 1.00     2187
treatmentPC2                    0.10      0.35    -0.57     0.79 1.00     1728
date2024M01M12:treatmentMC1    -0.27      0.33    -0.92     0.38 1.00     2988
date2024M01M23:treatmentMC1    -0.32      0.33    -0.96     0.33 1.00     2830
date2024M02M07:treatmentMC1    -0.55      0.32    -1.20     0.07 1.00     2961
date2024M02M21:treatmentMC1    -0.73      0.33    -1.39    -0.09 1.00     2731
date2024M03M06:treatmentMC1    -0.81      0.33    -1.47    -0.15 1.00     2874
date2024M03M22:treatmentMC1    -0.80      0.33    -1.44    -0.15 1.00     2825
date2024M04M04:treatmentMC1    -0.92      0.35    -1.60    -0.23 1.00     2874
date2024M04M16:treatmentMC1    -1.01      0.33    -1.66    -0.37 1.00     2227
date2025M02M08:treatmentMC1     0.01      0.34    -0.67     0.66 1.00     2500
date2025M02M15:treatmentMC1    -0.07      0.33    -0.72     0.57 1.00     2770
date2025M03M04:treatmentMC1    -0.10      0.34    -0.78     0.56 1.00     3144
date2025M03M16:treatmentMC1    -0.47      0.32    -1.12     0.16 1.00     2787
date2025M03M27:treatmentMC1    -0.50      0.34    -1.19     0.16 1.00     2750
date2025M04M11:treatmentMC1    -0.40      0.36    -1.11     0.29 1.00     3045
date2024M01M12:treatmentMC2    -0.04      0.33    -0.69     0.62 1.00     2875
date2024M01M23:treatmentMC2     0.06      0.33    -0.59     0.70 1.00     2775
date2024M02M07:treatmentMC2    -0.19      0.33    -0.84     0.45 1.00     2798
date2024M02M21:treatmentMC2    -0.14      0.34    -0.82     0.51 1.00     2710
date2024M03M06:treatmentMC2    -0.10      0.34    -0.76     0.57 1.00     2825
date2024M03M22:treatmentMC2    -0.11      0.33    -0.77     0.53 1.00     2511
date2024M04M04:treatmentMC2    -0.19      0.34    -0.86     0.48 1.00     2886
date2024M04M16:treatmentMC2    -0.37      0.34    -1.04     0.29 1.00     2516
date2025M02M08:treatmentMC2     0.19      0.35    -0.49     0.86 1.00     2463
date2025M02M15:treatmentMC2     0.66      0.34    -0.01     1.31 1.00     2887
date2025M03M04:treatmentMC2     0.55      0.35    -0.15     1.23 1.00     3094
date2025M03M16:treatmentMC2     0.21      0.33    -0.44     0.85 1.00     2700
date2025M03M27:treatmentMC2     0.12      0.34    -0.56     0.79 1.00     3062
date2025M04M11:treatmentMC2     0.12      0.39    -0.64     0.88 1.00     3304
date2024M01M12:treatmentPC0    -0.31      0.33    -0.95     0.35 1.00     2992
date2024M01M23:treatmentPC0    -0.56      0.32    -1.19     0.07 1.00     2745
date2024M02M07:treatmentPC0    -0.40      0.33    -1.05     0.23 1.00     2922
date2024M02M21:treatmentPC0    -0.37      0.33    -1.02     0.28 1.00     3116
date2024M03M06:treatmentPC0    -0.36      0.33    -1.00     0.30 1.00     2483
date2024M03M22:treatmentPC0    -0.54      0.33    -1.19     0.11 1.00     2706
date2024M04M04:treatmentPC0    -0.59      0.34    -1.25     0.09 1.00     2953
date2024M04M16:treatmentPC0    -0.51      0.33    -1.17     0.13 1.00     3257
date2025M02M08:treatmentPC0    -0.59      0.39    -1.35     0.17 1.00     3075
date2025M02M15:treatmentPC0    -0.55      0.33    -1.21     0.10 1.00     2969
date2025M03M04:treatmentPC0    -0.46      0.36    -1.17     0.26 1.00     3014
date2025M03M16:treatmentPC0    -0.81      0.34    -1.47    -0.14 1.00     3140
date2025M03M27:treatmentPC0    -0.67      0.38    -1.39     0.08 1.00     3335
date2025M04M11:treatmentPC0    -0.53      0.37    -1.25     0.19 1.00     2247
date2024M01M12:treatmentPC1    -0.46      0.33    -1.11     0.19 1.00     2995
date2024M01M23:treatmentPC1    -0.60      0.33    -1.23     0.05 1.00     2892
date2024M02M07:treatmentPC1    -0.91      0.32    -1.54    -0.29 1.00     2752
date2024M02M21:treatmentPC1    -0.84      0.33    -1.48    -0.21 1.00     2476
date2024M03M06:treatmentPC1    -0.73      0.33    -1.37    -0.08 1.00     2880
date2024M03M22:treatmentPC1    -0.82      0.33    -1.46    -0.18 1.00     2772
date2024M04M04:treatmentPC1    -0.95      0.34    -1.61    -0.28 1.00     2979
date2024M04M16:treatmentPC1    -1.03      0.34    -1.70    -0.37 1.00     2945
date2025M02M08:treatmentPC1    -1.01      0.33    -1.66    -0.36 1.00     2785
date2025M02M15:treatmentPC1    -1.05      0.34    -1.72    -0.39 1.00     2882
date2025M03M04:treatmentPC1    -1.12      0.34    -1.80    -0.45 1.00     3182
date2025M03M16:treatmentPC1    -1.05      0.34    -1.71    -0.39 1.00     3087
date2025M03M27:treatmentPC1    -1.04      0.35    -1.72    -0.34 1.00     3212
date2025M04M11:treatmentPC1    -1.23      0.40    -2.00    -0.44 1.00     3125
date2024M01M12:treatmentPC2    -0.29      0.36    -1.00     0.43 1.00     2778
date2024M01M23:treatmentPC2    -0.55      0.36    -1.25     0.14 1.00     2333
date2024M02M07:treatmentPC2    -0.78      0.36    -1.47    -0.07 1.00     2159
date2024M02M21:treatmentPC2    -0.61      0.36    -1.31     0.11 1.00     2491
date2024M03M06:treatmentPC2    -0.58      0.37    -1.30     0.13 1.00     1856
date2024M03M22:treatmentPC2    -0.55      0.36    -1.26     0.15 1.00     2027
date2024M04M04:treatmentPC2    -0.62      0.37    -1.34     0.12 1.00     2706
date2024M04M16:treatmentPC2    -0.68      0.36    -1.39     0.02 1.00     3141
date2025M02M08:treatmentPC2    -0.65      0.37    -1.37     0.07 1.00     2916
date2025M02M15:treatmentPC2    -0.36      0.36    -1.07     0.35 1.00     2638
date2025M03M04:treatmentPC2    -0.22      0.38    -0.96     0.52 1.00     2541
date2025M03M16:treatmentPC2    -0.40      0.36    -1.11     0.30 1.00     2830
date2025M03M27:treatmentPC2    -0.38      0.38    -1.14     0.36 1.00     2334
date2025M04M11:treatmentPC2    -0.56      0.43    -1.38     0.30 1.00     2570
                            Tail_ESS
Intercept                        489
date2024M01M12                  4611
date2024M01M23                  4161
date2024M02M07                  4299
date2024M02M21                  5113
date2024M03M06                  5267
date2024M03M22                  5165
date2024M04M04                  5257
date2024M04M16                  5252
date2025M02M08                  4562
date2025M02M15                  5271
date2025M03M04                  5281
date2025M03M16                  5087
date2025M03M27                  5755
date2025M04M11                  5094
treatmentMC1                    5223
treatmentMC2                    6142
treatmentPC0                    6150
treatmentPC1                    5902
treatmentPC2                    1876
date2024M01M12:treatmentMC1     7291
date2024M01M23:treatmentMC1     6716
date2024M02M07:treatmentMC1     6611
date2024M02M21:treatmentMC1     5988
date2024M03M06:treatmentMC1     6631
date2024M03M22:treatmentMC1     6946
date2024M04M04:treatmentMC1     6233
date2024M04M16:treatmentMC1     5203
date2025M02M08:treatmentMC1     6029
date2025M02M15:treatmentMC1     5923
date2025M03M04:treatmentMC1     6228
date2025M03M16:treatmentMC1     5853
date2025M03M27:treatmentMC1     7492
date2025M04M11:treatmentMC1     7125
date2024M01M12:treatmentMC2     7081
date2024M01M23:treatmentMC2     6424
date2024M02M07:treatmentMC2     6043
date2024M02M21:treatmentMC2     6387
date2024M03M06:treatmentMC2     7069
date2024M03M22:treatmentMC2     5811
date2024M04M04:treatmentMC2     6399
date2024M04M16:treatmentMC2     7370
date2025M02M08:treatmentMC2     6477
date2025M02M15:treatmentMC2     6186
date2025M03M04:treatmentMC2     6360
date2025M03M16:treatmentMC2     6414
date2025M03M27:treatmentMC2     7132
date2025M04M11:treatmentMC2     6804
date2024M01M12:treatmentPC0     6866
date2024M01M23:treatmentPC0     7129
date2024M02M07:treatmentPC0     6789
date2024M02M21:treatmentPC0     6619
date2024M03M06:treatmentPC0     5943
date2024M03M22:treatmentPC0     7353
date2024M04M04:treatmentPC0     6545
date2024M04M16:treatmentPC0     7638
date2025M02M08:treatmentPC0     8050
date2025M02M15:treatmentPC0     6781
date2025M03M04:treatmentPC0     7508
date2025M03M16:treatmentPC0     6850
date2025M03M27:treatmentPC0     7941
date2025M04M11:treatmentPC0     5367
date2024M01M12:treatmentPC1     6645
date2024M01M23:treatmentPC1     6957
date2024M02M07:treatmentPC1     7158
date2024M02M21:treatmentPC1     6270
date2024M03M06:treatmentPC1     7322
date2024M03M22:treatmentPC1     6564
date2024M04M04:treatmentPC1     6755
date2024M04M16:treatmentPC1     7499
date2025M02M08:treatmentPC1     6374
date2025M02M15:treatmentPC1     7465
date2025M03M04:treatmentPC1     7572
date2025M03M16:treatmentPC1     7202
date2025M03M27:treatmentPC1     7874
date2025M04M11:treatmentPC1     8764
date2024M01M12:treatmentPC2     6428
date2024M01M23:treatmentPC2     6396
date2024M02M07:treatmentPC2     6143
date2024M02M21:treatmentPC2     6478
date2024M03M06:treatmentPC2     3113
date2024M03M22:treatmentPC2     4921
date2024M04M04:treatmentPC2     7305
date2024M04M16:treatmentPC2     6883
date2025M02M08:treatmentPC2     7278
date2025M02M15:treatmentPC2     7244
date2025M03M04:treatmentPC2     6715
date2025M03M16:treatmentPC2     6105
date2025M03M27:treatmentPC2     3233
date2025M04M11:treatmentPC2     7080

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     5.78      0.35     5.11     6.48 1.00    13355    11217

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fitted(fit_glmer_gamma_mg), residuals(fit_glmer_gamma_mg),
     xlab = "Fitted Values",
     ylab = "Residuals",
     main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)

emmip(fit_glmer_gamma_mg, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data %>% filter(study_year == 1),
             aes(x = date, y = mg, color = treatment),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  theme_bw()

Warning: Removed 32 rows containing missing values or values outside the scale range
(`geom_point()`).

emmip(fit_mg, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data %>% filter(study_year == 1),
             aes(x = date, y = mg, color = treatment),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  theme_bw()

Warning: Removed 32 rows containing missing values or values outside the scale range
(`geom_point()`).

P

summary(fit_p)

Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.

Warning: There were 3210 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: p | cens(p_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id) 
   Data: lc_model_data (Number of observations: 626) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~site_block_id (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.15      0.12     0.01     0.46 1.01      460      632

~trt_site_block_id (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.59      0.10     0.42     0.83 1.00      534      947

Regression Coefficients:
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                      -1.23      0.59    -2.34    -0.13 1.09       38
date2024M01M12                 -0.51      0.71    -1.89     0.84 1.10       36
date2024M01M23                 -1.28      0.71    -2.68     0.03 1.09       46
date2024M02M07                 -2.05      0.70    -3.39    -0.74 1.08       44
date2024M02M21                 -2.69      0.73    -4.07    -1.29 1.11       32
date2024M03M06                 -2.22      0.73    -3.69    -0.83 1.06       53
date2024M03M22                 -1.31      0.73    -2.73     0.20 1.08       41
date2024M04M04                 -1.37      0.74    -2.77     0.11 1.08       50
date2024M04M16                 -2.19      0.71    -3.51    -0.79 1.07       51
date2025M02M08                 -0.13      0.70    -1.49     1.20 1.07       55
date2025M02M15                 -1.31      0.75    -2.80     0.13 1.06       59
date2025M03M04               -284.60    205.13  -678.11   -11.41 1.02      117
date2025M03M16                 -1.42      0.76    -2.95     0.06 1.09       36
date2025M03M27                 -0.92      0.77    -2.36     0.61 1.06       58
date2025M04M11                 -0.07      0.87    -1.74     1.70 1.07       49
treatmentMC1                   -0.15      0.75    -1.59     1.31 1.09       37
treatmentMC2                   -0.36      0.74    -1.80     1.13 1.09       37
treatmentPC0                    0.54      0.80    -0.93     2.21 1.06       54
treatmentPC1                    0.14      0.76    -1.32     1.61 1.06       57
treatmentPC2                    0.18      0.93    -1.43     2.20 1.05       62
date2024M01M12:treatmentMC1    -0.25      0.95    -2.12     1.65 1.09       39
date2024M01M23:treatmentMC1     0.07      0.95    -1.78     1.92 1.08       51
date2024M02M07:treatmentMC1    -0.36      0.94    -2.20     1.54 1.07       48
date2024M02M21:treatmentMC1     0.37      0.96    -1.53     2.22 1.10       33
date2024M03M06:treatmentMC1    -0.49      0.95    -2.31     1.46 1.07       54
date2024M03M22:treatmentMC1    -1.74      0.98    -3.65     0.24 1.08       36
date2024M04M04:treatmentMC1    -1.41      0.96    -3.30     0.48 1.08       50
date2024M04M16:treatmentMC1    -0.62      0.94    -2.53     1.14 1.06       53
date2025M02M08:treatmentMC1     0.32      0.98    -1.58     2.16 1.06       68
date2025M02M15:treatmentMC1    -0.38      1.00    -2.26     1.58 1.05       65
date2025M03M04:treatmentMC1   282.26    205.12     8.99   675.95 1.02      117
date2025M03M16:treatmentMC1  -233.74    186.65  -650.85    -8.95 1.00      490
date2025M03M27:treatmentMC1  -235.65    186.00  -654.17    -9.65 1.00      496
date2025M04M11:treatmentMC1    -0.44      1.14    -2.65     1.75 1.06       55
date2024M01M12:treatmentMC2    -0.03      0.93    -1.81     1.80 1.07       56
date2024M01M23:treatmentMC2    -0.37      0.92    -2.15     1.41 1.07       59
date2024M02M07:treatmentMC2    -0.57      0.92    -2.40     1.21 1.07       58
date2024M02M21:treatmentMC2     0.19      0.91    -1.63     1.97 1.08       43
date2024M03M06:treatmentMC2    -0.53      0.96    -2.39     1.32 1.05       72
date2024M03M22:treatmentMC2    -1.96      0.93    -3.87    -0.18 1.07       49
date2024M04M04:treatmentMC2    -1.18      0.93    -3.01     0.63 1.07       58
date2024M04M16:treatmentMC2     0.64      0.96    -1.20     2.52 1.06       60
date2025M02M08:treatmentMC2    -0.87      0.94    -2.70     0.94 1.05       75
date2025M02M15:treatmentMC2     0.14      0.97    -1.73     2.04 1.05       67
date2025M03M04:treatmentMC2    31.68    275.55  -512.93   561.27 1.01      201
date2025M03M16:treatmentMC2  -230.00    181.27  -644.10    -8.83 1.00      572
date2025M03M27:treatmentMC2    -0.13      1.01    -2.16     1.79 1.05       68
date2025M04M11:treatmentMC2    -0.96      1.20    -3.26     1.46 1.04       85
date2024M01M12:treatmentPC0    -1.08      0.97    -3.03     0.77 1.05       69
date2024M01M23:treatmentPC0    -1.33      0.97    -3.25     0.50 1.05       65
date2024M02M07:treatmentPC0    -1.33      0.99    -3.21     0.56 1.04       76
date2024M02M21:treatmentPC0    -0.00      1.01    -1.98     1.87 1.05       75
date2024M03M06:treatmentPC0    -0.97      1.02    -2.98     1.04 1.04       67
date2024M03M22:treatmentPC0    -2.23      1.00    -4.25    -0.36 1.05       74
date2024M04M04:treatmentPC0    -1.29      1.00    -3.26     0.59 1.05       77
date2024M04M16:treatmentPC0    -0.25      0.97    -2.15     1.58 1.04       76
date2025M02M08:treatmentPC0    -0.11      1.20    -2.42     2.23 1.03       96
date2025M02M15:treatmentPC0     0.61      1.04    -1.44     2.56 1.04       67
date2025M03M04:treatmentPC0   283.09    205.16     9.75   676.78 1.02      117
date2025M03M16:treatmentPC0  -258.26    191.44  -658.91   -10.67 1.00      775
date2025M03M27:treatmentPC0    -0.10      1.15    -2.34     2.17 1.03      115
date2025M04M11:treatmentPC0    -0.06      1.22    -2.49     2.27 1.04       81
date2024M01M12:treatmentPC1    -0.27      0.93    -2.05     1.54 1.06       70
date2024M01M23:treatmentPC1    -0.51      0.94    -2.33     1.39 1.06       66
date2024M02M07:treatmentPC1     0.01      0.90    -1.75     1.72 1.05       73
date2024M02M21:treatmentPC1     0.27      0.93    -1.54     2.10 1.07       60
date2024M03M06:treatmentPC1    -0.37      0.94    -2.14     1.47 1.05       68
date2024M03M22:treatmentPC1    -1.56      0.94    -3.41     0.28 1.05       82
date2024M04M04:treatmentPC1    -0.95      0.96    -2.82     0.94 1.05       86
date2024M04M16:treatmentPC1     1.19      0.97    -0.73     3.07 1.04       90
date2025M02M08:treatmentPC1     1.25      0.97    -0.63     3.17 1.04       89
date2025M02M15:treatmentPC1     1.00      0.98    -0.93     2.93 1.03       86
date2025M03M04:treatmentPC1   283.12    205.13     9.83   676.79 1.02      117
date2025M03M16:treatmentPC1     0.64      0.99    -1.30     2.60 1.06       79
date2025M03M27:treatmentPC1     0.74      1.04    -1.32     2.75 1.04       92
date2025M04M11:treatmentPC1     1.19      1.23    -1.19     3.61 1.03       99
date2024M01M12:treatmentPC2    -0.25      1.07    -2.43     1.71 1.06       55
date2024M01M23:treatmentPC2     0.11      1.06    -2.12     2.04 1.06       53
date2024M02M07:treatmentPC2     0.33      1.06    -1.88     2.22 1.06       58
date2024M02M21:treatmentPC2    -0.03      1.06    -2.13     1.93 1.06       47
date2024M03M06:treatmentPC2    -0.33      1.06    -2.53     1.64 1.04       69
date2024M03M22:treatmentPC2    -0.63      1.08    -2.88     1.35 1.05       58
date2024M04M04:treatmentPC2     0.25      1.07    -1.98     2.27 1.05       71
date2024M04M16:treatmentPC2     0.65      1.08    -1.56     2.63 1.05       71
date2025M02M08:treatmentPC2     0.71      1.08    -1.47     2.79 1.05       69
date2025M02M15:treatmentPC2     1.23      1.07    -0.91     3.24 1.04       78
date2025M03M04:treatmentPC2   283.98    205.09    10.88   677.53 1.02      117
date2025M03M16:treatmentPC2     1.42      1.09    -0.76     3.51 1.05       56
date2025M03M27:treatmentPC2     1.38      1.16    -1.04     3.58 1.04       74
date2025M04M11:treatmentPC2     1.30      1.33    -1.33     3.90 1.04       74
                            Tail_ESS
Intercept                        152
date2024M01M12                   166
date2024M01M23                   253
date2024M02M07                   192
date2024M02M21                   351
date2024M03M06                    96
date2024M03M22                   327
date2024M04M04                   243
date2024M04M16                   184
date2025M02M08                   312
date2025M02M15                   106
date2025M03M04                   366
date2025M03M16                   100
date2025M03M27                   558
date2025M04M11                   135
treatmentMC1                     231
treatmentMC2                     184
treatmentPC0                     194
treatmentPC1                     299
treatmentPC2                     154
date2024M01M12:treatmentMC1      202
date2024M01M23:treatmentMC1      282
date2024M02M07:treatmentMC1      223
date2024M02M21:treatmentMC1      317
date2024M03M06:treatmentMC1      129
date2024M03M22:treatmentMC1      338
date2024M04M04:treatmentMC1      264
date2024M04M16:treatmentMC1      240
date2025M02M08:treatmentMC1      385
date2025M02M15:treatmentMC1      211
date2025M03M04:treatmentMC1      361
date2025M03M16:treatmentMC1     1068
date2025M03M27:treatmentMC1      956
date2025M04M11:treatmentMC1      157
date2024M01M12:treatmentMC2      346
date2024M01M23:treatmentMC2      303
date2024M02M07:treatmentMC2      248
date2024M02M21:treatmentMC2      216
date2024M03M06:treatmentMC2      176
date2024M03M22:treatmentMC2      316
date2024M04M04:treatmentMC2      511
date2024M04M16:treatmentMC2      595
date2025M02M08:treatmentMC2      377
date2025M02M15:treatmentMC2      210
date2025M03M04:treatmentMC2      670
date2025M03M16:treatmentMC2      993
date2025M03M27:treatmentMC2      590
date2025M04M11:treatmentMC2      579
date2024M01M12:treatmentPC0      310
date2024M01M23:treatmentPC0      232
date2024M02M07:treatmentPC0      280
date2024M02M21:treatmentPC0      407
date2024M03M06:treatmentPC0      186
date2024M03M22:treatmentPC0      281
date2024M04M04:treatmentPC0      314
date2024M04M16:treatmentPC0      183
date2025M02M08:treatmentPC0      471
date2025M02M15:treatmentPC0      254
date2025M03M04:treatmentPC0      366
date2025M03M16:treatmentPC0     1263
date2025M03M27:treatmentPC0      654
date2025M04M11:treatmentPC0      478
date2024M01M12:treatmentPC1      413
date2024M01M23:treatmentPC1      323
date2024M02M07:treatmentPC1      462
date2024M02M21:treatmentPC1      374
date2024M03M06:treatmentPC1      213
date2024M03M22:treatmentPC1      400
date2024M04M04:treatmentPC1      388
date2024M04M16:treatmentPC1      405
date2025M02M08:treatmentPC1      365
date2025M02M15:treatmentPC1      265
date2025M03M04:treatmentPC1      348
date2025M03M16:treatmentPC1      560
date2025M03M27:treatmentPC1      679
date2025M04M11:treatmentPC1      559
date2024M01M12:treatmentPC2      177
date2024M01M23:treatmentPC2      203
date2024M02M07:treatmentPC2      168
date2024M02M21:treatmentPC2      186
date2024M03M06:treatmentPC2      231
date2024M03M22:treatmentPC2      226
date2024M04M04:treatmentPC2      315
date2024M04M16:treatmentPC2      260
date2025M02M08:treatmentPC2      290
date2025M02M15:treatmentPC2      205
date2025M03M04:treatmentPC2      363
date2025M03M16:treatmentPC2      246
date2025M03M27:treatmentPC2      212
date2025M04M11:treatmentPC2      612

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.80      0.06     0.69     0.91 1.00      863     1421

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

# neff_ratio(fit_p)
# prior_summary(fit_p)
# plot(fit, type = "trace", nvariables = 1)

pp_check(fit_p) + xlim(0, 1)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 189 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 18 rows containing non-finite outside the scale range
(`stat_density()`).

pp_check(fit_p, type = "dens_overlay", ndraws = 100)  + xlim(0, 1)

Warning: Censored responses are not shown in 'pp_check'.

Warning: Removed 1818 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 18 rows containing non-finite outside the scale range
(`stat_density()`).

# loo(fit_p)

conditional_effects(fit_p, effects = "treatment")

conditional_effects(fit_p, effects = "date:treatment")

emm_treatments <- emmeans(fit_p, ~ treatment, type="response")

NOTE: Results may be misleading due to involvement in interactions

emm_treatments

 treatment response lower.HPD upper.HPD
 MC0       1.00e-08    0.0000  2.37e-02
 MC1       0.00e+00    0.0000  8.72e-05
 MC2       0.00e+00    0.0000  4.81e-05
 PC0       5.00e-08    0.0000  2.28e-02
 PC1       1.12e-01    0.0620  1.76e-01
 PC2       1.57e-01    0.0831  2.54e-01

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

pairs(emm_treatments)

 contrast     ratio lower.HPD upper.HPD
 MC0 / MC1 146514.1     0.000  2.59e+22
 MC0 / MC2 463834.4     0.000  1.79e+22
 MC0 / PC0      0.3     0.000  3.82e+12
 MC0 / PC1      0.0     0.000  0.00e+00
 MC0 / PC2      0.0     0.000  0.00e+00
 MC1 / MC2      3.0     0.000  4.38e+17
 MC1 / PC0      0.0     0.000  2.34e+08
 MC1 / PC1      0.0     0.000  0.00e+00
 MC1 / PC2      0.0     0.000  0.00e+00
 MC2 / PC0      0.0     0.000  7.10e+07
 MC2 / PC1      0.0     0.000  0.00e+00
 MC2 / PC2      0.0     0.000  0.00e+00
 PC0 / PC1      0.0     0.000  0.00e+00
 PC0 / PC2      0.0     0.000  0.00e+00
 PC1 / PC2      0.7     0.306  1.00e+00

Results are averaged over the levels of: date 
Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

plot(emm_treatments) + theme_bw()

emmip(fit_p, ~ date | treatment, type = "response", CIs = T) + 
  geom_point(data = lc_model_data,
             aes(x = date, y = p, color = as.factor(p_censored)),
             alpha = 0.7,
             position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
  # ylim(0, 4) +
  theme_bw()

Warning: Removed 94 rows containing missing values or values outside the scale range
(`geom_point()`).

Zero-Inflated

Nitrate

fit_no3 <- glmmTMB(
    no3 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
    zi = ~ treatment,
    family = glmmTMB::lognormal(link = "log"),
    data = df_model
)

Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

S(fit_no3)

 Family: lognormal  ( log )
Formula:          
no3 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation:       ~treatment
Data: df_model

      AIC       BIC    logLik -2*log(L)  df.resid 
       NA        NA        NA        NA       343 

Random effects:

Conditional model:
 Groups            Name        Variance Std.Dev.
 site_block_id     (Intercept) 0.18431  0.4293  
 trt_site_block_id (Intercept) 0.07486  0.2736  
Number of obs: 406, groups:  site_block_id, 8; trt_site_block_id, 48

Dispersion parameter for lognormal family (): 43.2 

Conditional model:
                            Estimate Std. Error z value Pr(>|z|)
(Intercept)                  4.94296        NaN     NaN      NaN
date2024-01-12              -0.46589        NaN     NaN      NaN
date2024-01-23              -2.24045        NaN     NaN      NaN
date2024-02-07              -1.16452        NaN     NaN      NaN
date2024-02-21              -0.93779        NaN     NaN      NaN
date2024-03-06              -0.58716        NaN     NaN      NaN
date2024-03-22              -2.06985        NaN     NaN      NaN
date2024-04-04              -0.85750        NaN     NaN      NaN
date2024-04-16              -0.35219        NaN     NaN      NaN
treatmentMC1                -0.95229        NaN     NaN      NaN
treatmentMC2                -0.18415        NaN     NaN      NaN
treatmentPC0                -0.59872        NaN     NaN      NaN
treatmentPC1                -0.43191        NaN     NaN      NaN
treatmentPC2                -0.26465        NaN     NaN      NaN
date2024-01-12:treatmentMC1  0.35628        NaN     NaN      NaN
date2024-01-23:treatmentMC1  2.13208        NaN     NaN      NaN
date2024-02-07:treatmentMC1 -0.03654        NaN     NaN      NaN
date2024-02-21:treatmentMC1  0.00000        NaN     NaN      NaN
date2024-03-06:treatmentMC1  0.00000        NaN     NaN      NaN
date2024-03-22:treatmentMC1  0.71224        NaN     NaN      NaN
date2024-04-04:treatmentMC1  0.00000        NaN     NaN      NaN
date2024-04-16:treatmentMC1  0.00000        NaN     NaN      NaN
date2024-01-12:treatmentMC2 -0.51072        NaN     NaN      NaN
date2024-01-23:treatmentMC2  0.20048        NaN     NaN      NaN
date2024-02-07:treatmentMC2  0.48826        NaN     NaN      NaN
date2024-02-21:treatmentMC2 -0.29661        NaN     NaN      NaN
date2024-03-06:treatmentMC2 -1.29413        NaN     NaN      NaN
date2024-03-22:treatmentMC2  0.00000        NaN     NaN      NaN
date2024-04-04:treatmentMC2  0.00000        NaN     NaN      NaN
date2024-04-16:treatmentMC2 -1.03794        NaN     NaN      NaN
date2024-01-12:treatmentPC0 -0.12570        NaN     NaN      NaN
date2024-01-23:treatmentPC0  0.84651        NaN     NaN      NaN
date2024-02-07:treatmentPC0  0.00000        NaN     NaN      NaN
date2024-02-21:treatmentPC0 -1.56917        NaN     NaN      NaN
date2024-03-06:treatmentPC0  0.00000        NaN     NaN      NaN
date2024-03-22:treatmentPC0  0.00000        NaN     NaN      NaN
date2024-04-04:treatmentPC0  0.00000        NaN     NaN      NaN
date2024-04-16:treatmentPC0  0.00000        NaN     NaN      NaN
date2024-01-12:treatmentPC1 -0.32762        NaN     NaN      NaN
date2024-01-23:treatmentPC1  0.64563        NaN     NaN      NaN
date2024-02-07:treatmentPC1  1.88958        NaN     NaN      NaN
date2024-02-21:treatmentPC1  0.92799        NaN     NaN      NaN
date2024-03-06:treatmentPC1  1.48469        NaN     NaN      NaN
date2024-03-22:treatmentPC1  1.94375        NaN     NaN      NaN
date2024-04-04:treatmentPC1  0.00000        NaN     NaN      NaN
date2024-04-16:treatmentPC1  0.68575        NaN     NaN      NaN
date2024-01-12:treatmentPC2 -0.11587        NaN     NaN      NaN
date2024-01-23:treatmentPC2  0.56902        NaN     NaN      NaN
date2024-02-07:treatmentPC2  0.00000        NaN     NaN      NaN
date2024-02-21:treatmentPC2  0.00000        NaN     NaN      NaN
date2024-03-06:treatmentPC2 -0.77773        NaN     NaN      NaN
date2024-03-22:treatmentPC2 -0.03610        NaN     NaN      NaN
date2024-04-04:treatmentPC2 -0.85750        NaN     NaN      NaN
date2024-04-16:treatmentPC2  0.00000        NaN     NaN      NaN

Zero-inflation model:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    0.6702        NaN     NaN      NaN
treatmentMC1   0.1565        NaN     NaN      NaN
treatmentMC2   0.2974        NaN     NaN      NaN
treatmentPC0   0.3366        NaN     NaN      NaN
treatmentPC1   0.3313        NaN     NaN      NaN
treatmentPC2   0.2141        NaN     NaN      NaN

simulationOutput <- simulateResiduals(fittedModel = fit_no3, n = 500)
plot(simulationOutput)

testUniformity(simulationOutput)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.035847, p-value = 0.6739
alternative hypothesis: two-sided

testOutliers(simulationOutput)

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 2, observations = 406, p-value = 0.6793
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
 0.0005971323 0.0176806666
sample estimates:
frequency of outliers (expected: 0.00399201596806387 ) 
                                           0.004926108 

testDispersion(simulationOutput)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.84685, p-value = 1
alternative hypothesis: two.sided

testZeroInflation(simulationOutput)

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 0.99872, p-value = 1
alternative hypothesis: two.sided

library(ggeffects)

Warning: package 'ggeffects' was built under R version 4.4.3

library(RColorBrewer)
pred_effects <- ggpredict(fit_no3, terms = c("date", "treatment"))

You are calculating adjusted predictions on the population-level (i.e.
  `type = "fixed"`) for a *generalized* linear mixed model.
  This may produce biased estimates due to Jensen's inequality. Consider
  setting `bias_correction = TRUE` to correct for this bias.
  See also the documentation of the `bias_correction` argument.

plot(pred_effects, colors = "flat") +
  theme_bw() +
  labs(
    x = "Date",
    y = "Nitrate",
    title = "Predicted Values of Nitrate",
    color = "Treatment"
  )

print(pred_effects)

# Predicted values of no3

treatment: MC0

date       | Predicted
----------------------
2023-12-27 |    140.18
2024-01-12 |     87.98
2024-01-23 |     14.92
2024-02-07 |     43.75
2024-02-21 |     54.88
2024-03-06 |     77.93
2024-03-22 |     17.69
2024-04-04 |     59.47
2024-04-16 |     98.57

treatment: MC1

date       | Predicted
----------------------
2023-12-27 |     54.09
2024-01-12 |     48.48
2024-01-23 |     48.54
2024-02-07 |     16.27
2024-02-21 |     21.18
2024-03-06 |     30.07
2024-03-22 |     13.92
2024-04-04 |     22.95
2024-04-16 |     38.03

treatment: MC2

date       | Predicted
----------------------
2023-12-27 |    116.61
2024-01-12 |     43.91
2024-01-23 |     15.16
2024-02-07 |     59.30
2024-02-21 |     33.93
2024-03-06 |     17.77
2024-03-22 |     14.72
2024-04-04 |     49.47
2024-04-16 |     29.04

treatment: PC0

date       | Predicted
----------------------
2023-12-27 |     77.03
2024-01-12 |     42.63
2024-01-23 |     19.11
2024-02-07 |     24.04
2024-02-21 |      6.28
2024-03-06 |     42.82
2024-03-22 |      9.72
2024-04-04 |     32.68
2024-04-16 |     54.17

treatment: PC1

date       | Predicted
----------------------
2023-12-27 |     91.02
2024-01-12 |     41.16
2024-01-23 |     18.47
2024-02-07 |    187.94
2024-02-21 |     90.13
2024-03-06 |    223.32
2024-03-22 |     80.23
2024-04-04 |     38.61
2024-04-16 |    127.05

treatment: PC2

date       | Predicted
----------------------
2023-12-27 |    107.59
2024-01-12 |     60.13
2024-01-23 |     20.22
2024-02-07 |     33.58
2024-02-21 |     42.12
2024-03-06 |     27.48
2024-03-22 |     13.10
2024-04-04 |     19.36
2024-04-16 |     75.65

Adjusted for:
*     site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)

NH4

fit_nh4 <- glmmTMB(
    nh4 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
    zi = ~ treatment,
    family = glmmTMB::lognormal(link = "log"),
    data = df_model
)

Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

S(fit_nh4)

 Family: lognormal  ( log )
Formula:          
nh4 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation:       ~treatment
Data: df_model

      AIC       BIC    logLik -2*log(L)  df.resid 
       NA        NA        NA        NA       349 

Random effects:

Conditional model:
 Groups            Name        Variance  Std.Dev. 
 site_block_id     (Intercept) 4.380e-11 6.618e-06
 trt_site_block_id (Intercept) 3.594e-10 1.896e-05
Number of obs: 412, groups:  site_block_id, 8; trt_site_block_id, 48

Dispersion parameter for lognormal family (): 1.13 

Conditional model:
                             Estimate Std. Error z value Pr(>|z|)
(Intercept)                  0.543717        NaN     NaN      NaN
date2024-01-12              -0.606122        NaN     NaN      NaN
date2024-01-23              -0.068191        NaN     NaN      NaN
date2024-02-07              -0.119794        NaN     NaN      NaN
date2024-02-21              -0.365954        NaN     NaN      NaN
date2024-03-06              -0.386029        NaN     NaN      NaN
date2024-03-22              -0.702810        NaN     NaN      NaN
date2024-04-04              -0.346473        NaN     NaN      NaN
date2024-04-16              -0.444480        NaN     NaN      NaN
treatmentMC1                -0.186478        NaN     NaN      NaN
treatmentMC2                -0.300228        NaN     NaN      NaN
treatmentPC0                -0.438777        NaN     NaN      NaN
treatmentPC1                -0.017852        NaN     NaN      NaN
treatmentPC2                -0.568562        NaN     NaN      NaN
date2024-01-12:treatmentMC1 -0.009409        NaN     NaN      NaN
date2024-01-23:treatmentMC1 -0.475401        NaN     NaN      NaN
date2024-02-07:treatmentMC1  0.000000        NaN     NaN      NaN
date2024-02-21:treatmentMC1  0.000000        NaN     NaN      NaN
date2024-03-06:treatmentMC1  0.000000        NaN     NaN      NaN
date2024-03-22:treatmentMC1  0.000000        NaN     NaN      NaN
date2024-04-04:treatmentMC1  0.331027        NaN     NaN      NaN
date2024-04-16:treatmentMC1  0.000000        NaN     NaN      NaN
date2024-01-12:treatmentMC2  0.362826        NaN     NaN      NaN
date2024-01-23:treatmentMC2 -0.689456        NaN     NaN      NaN
date2024-02-07:treatmentMC2  0.000000        NaN     NaN      NaN
date2024-02-21:treatmentMC2  0.000000        NaN     NaN      NaN
date2024-03-06:treatmentMC2  0.000000        NaN     NaN      NaN
date2024-03-22:treatmentMC2  0.000000        NaN     NaN      NaN
date2024-04-04:treatmentMC2  0.027561        NaN     NaN      NaN
date2024-04-16:treatmentMC2  0.000000        NaN     NaN      NaN
date2024-01-12:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-01-23:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-02-07:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-02-21:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-03-06:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-03-22:treatmentPC0  0.000000        NaN     NaN      NaN
date2024-04-04:treatmentPC0  0.026161        NaN     NaN      NaN
date2024-04-16:treatmentPC0  0.136907        NaN     NaN      NaN
date2024-01-12:treatmentPC1 -0.251759        NaN     NaN      NaN
date2024-01-23:treatmentPC1  0.889392        NaN     NaN      NaN
date2024-02-07:treatmentPC1  0.000000        NaN     NaN      NaN
date2024-02-21:treatmentPC1 -0.365954        NaN     NaN      NaN
date2024-03-06:treatmentPC1  0.000000        NaN     NaN      NaN
date2024-03-22:treatmentPC1  0.000000        NaN     NaN      NaN
date2024-04-04:treatmentPC1  0.112345        NaN     NaN      NaN
date2024-04-16:treatmentPC1  0.547386        NaN     NaN      NaN
date2024-01-12:treatmentPC2  0.275519        NaN     NaN      NaN
date2024-01-23:treatmentPC2  0.207274        NaN     NaN      NaN
date2024-02-07:treatmentPC2 -0.119794        NaN     NaN      NaN
date2024-02-21:treatmentPC2  0.000000        NaN     NaN      NaN
date2024-03-06:treatmentPC2  0.000000        NaN     NaN      NaN
date2024-03-22:treatmentPC2  0.000000        NaN     NaN      NaN
date2024-04-04:treatmentPC2  0.172446        NaN     NaN      NaN
date2024-04-16:treatmentPC2  0.000000        NaN     NaN      NaN

Zero-inflation model:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    1.2040        NaN     NaN      NaN
treatmentMC1   0.2741        NaN     NaN      NaN
treatmentMC2   0.5710        NaN     NaN      NaN
treatmentPC0   1.0245        NaN     NaN      NaN
treatmentPC1   0.2564        NaN     NaN      NaN
treatmentPC2   0.9445        NaN     NaN      NaN

simulationOutput <- simulateResiduals(fittedModel = fit_nh4, n = 500)
plot(simulationOutput)

testUniformity(simulationOutput)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.033539, p-value = 0.743
alternative hypothesis: two-sided

testOutliers(simulationOutput)

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 2, observations = 412, p-value = 0.6822
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
 0.0005884282 0.0174248139
sample estimates:
frequency of outliers (expected: 0.00399201596806387 ) 
                                           0.004854369 

testDispersion(simulationOutput)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.7099, p-value = 0.132
alternative hypothesis: two.sided

testZeroInflation(simulationOutput)

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 0.99968, p-value = 1
alternative hypothesis: two.sided

library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_nh4, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") + 
  theme_bw() +
  labs(
    x = "Date",
    y = "NH4",
    title = "Predicted Values of NH4",
    color = "Treatment"
  )

print(pred_effects)

# Predicted values of nh4

treatment: MC0

date       | Predicted
----------------------
2023-12-27 |      1.72
2024-01-12 |      0.94
2024-01-23 |      1.61
2024-02-07 |      1.53
2024-02-21 |      1.19
2024-03-06 |      1.17
2024-03-22 |      0.85
2024-04-04 |      1.22
2024-04-16 |      1.10

treatment: MC1

date       | Predicted
----------------------
2023-12-27 |      1.43
2024-01-12 |      0.77
2024-01-23 |      0.83
2024-02-07 |      1.27
2024-02-21 |      0.99
2024-03-06 |      0.97
2024-03-22 |      0.71
2024-04-04 |      1.41
2024-04-16 |      0.92

treatment: MC2

date       | Predicted
----------------------
2023-12-27 |      1.28
2024-01-12 |      1.00
2024-01-23 |      0.60
2024-02-07 |      1.13
2024-02-21 |      0.88
2024-03-06 |      0.87
2024-03-22 |      0.63
2024-04-04 |      0.93
2024-04-16 |      0.82

treatment: PC0

date       | Predicted
----------------------
2023-12-27 |      1.11
2024-01-12 |      0.61
2024-01-23 |      1.04
2024-02-07 |      0.99
2024-02-21 |      0.77
2024-03-06 |      0.75
2024-03-22 |      0.55
2024-04-04 |      0.81
2024-04-16 |      0.82

treatment: PC1

date       | Predicted
----------------------
2023-12-27 |      1.69
2024-01-12 |      0.72
2024-01-23 |      3.85
2024-02-07 |      1.50
2024-02-21 |      0.81
2024-03-06 |      1.15
2024-03-22 |      0.84
2024-04-04 |      1.34
2024-04-16 |      1.88

treatment: PC2

date       | Predicted
----------------------
2023-12-27 |      0.98
2024-01-12 |      0.70
2024-01-23 |      1.12
2024-02-07 |      0.77
2024-02-21 |      0.68
2024-03-06 |      0.66
2024-03-22 |      0.48
2024-04-04 |      0.82
2024-04-16 |      0.63

Adjusted for:
*     site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)

Mg

fit_mg <- glmmTMB(
    mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
    zi = ~ treatment,
    family = glmmTMB::lognormal(link = "log"),
    data = df_model
)

S(fit_mg)

 Family: lognormal  ( log )
Formula:          
mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation:      ~treatment
Data: df_model

      AIC       BIC    logLik -2*log(L)  df.resid 
   3712.5    3964.0   -1793.3    3586.5       337 

Random effects:

Conditional model:
 Groups            Name        Variance Std.Dev.
 site_block_id     (Intercept) 0.08437  0.2905  
 trt_site_block_id (Intercept) 0.06995  0.2645  
Number of obs: 400, groups:  site_block_id, 8; trt_site_block_id, 48

Dispersion parameter for lognormal family (): 20.5 

Conditional model:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                  3.89132    0.19848  19.605  < 2e-16 ***
date2024-01-12               0.19747    0.18513   1.067 0.286123    
date2024-01-23               0.31634    0.15943   1.984 0.047241 *  
date2024-02-07               0.12808    0.16715   0.766 0.443514    
date2024-02-21               0.21117    0.16575   1.274 0.202633    
date2024-03-06               0.32937    0.17429   1.890 0.058790 .  
date2024-03-22               0.48842    0.16919   2.887 0.003891 ** 
date2024-04-04               0.55742    0.17043   3.271 0.001073 ** 
date2024-04-16               0.41760    0.17030   2.452 0.014203 *  
treatmentMC1                 0.63274    0.21271   2.975 0.002934 ** 
treatmentMC2                 0.25160    0.22258   1.130 0.258313    
treatmentPC0                -0.66214    0.25529  -2.594 0.009497 ** 
treatmentPC1                -0.10683    0.25276  -0.423 0.672552    
treatmentPC2                -0.36107    0.28614  -1.262 0.206999    
date2024-01-12:treatmentMC1 -0.12862    0.21072  -0.610 0.541599    
date2024-01-23:treatmentMC1 -1.38917    0.28870  -4.812  1.5e-06 ***
date2024-02-07:treatmentMC1 -0.50278    0.20987  -2.396 0.016588 *  
date2024-02-21:treatmentMC1 -0.54368    0.20679  -2.629 0.008561 ** 
date2024-03-06:treatmentMC1 -0.64073    0.21213  -3.020 0.002524 ** 
date2024-03-22:treatmentMC1 -0.64668    0.20587  -3.141 0.001683 ** 
date2024-04-04:treatmentMC1 -0.66827    0.21603  -3.093 0.001978 ** 
date2024-04-16:treatmentMC1 -0.78931    0.21531  -3.666 0.000246 ***
date2024-01-12:treatmentMC2  0.10814    0.22548   0.480 0.631515    
date2024-01-23:treatmentMC2  0.10015    0.20230   0.495 0.620537    
date2024-02-07:treatmentMC2  0.02388    0.21279   0.112 0.910640    
date2024-02-21:treatmentMC2  0.10738    0.20869   0.515 0.606875    
date2024-03-06:treatmentMC2  0.13378    0.21369   0.626 0.531277    
date2024-03-22:treatmentMC2  0.13797    0.20790   0.664 0.506919    
date2024-04-04:treatmentMC2  0.05463    0.21106   0.259 0.795775    
date2024-04-16:treatmentMC2 -0.30114    0.22976  -1.311 0.189974    
date2024-01-12:treatmentPC0  0.25153    0.25958   0.969 0.332538    
date2024-01-23:treatmentPC0  0.28432    0.24273   1.171 0.241472    
date2024-02-07:treatmentPC0  0.36383    0.25736   1.414 0.157443    
date2024-02-21:treatmentPC0  0.41553    0.25202   1.649 0.099185 .  
date2024-03-06:treatmentPC0  0.39201    0.25491   1.538 0.124083    
date2024-03-22:treatmentPC0  0.28602    0.24940   1.147 0.251441    
date2024-04-04:treatmentPC0  0.28631    0.25015   1.145 0.252393    
date2024-04-16:treatmentPC0  0.35281    0.26061   1.354 0.175802    
date2024-01-12:treatmentPC1 -0.02233    0.26210  -0.085 0.932094    
date2024-01-23:treatmentPC1 -0.19236    0.25388  -0.758 0.448642    
date2024-02-07:treatmentPC1 -0.22277    0.25231  -0.883 0.377271    
date2024-02-21:treatmentPC1 -0.18247    0.25305  -0.721 0.470860    
date2024-03-06:treatmentPC1 -0.13884    0.25573  -0.543 0.587192    
date2024-03-22:treatmentPC1 -0.16120    0.25385  -0.635 0.525414    
date2024-04-04:treatmentPC1 -0.19426    0.26067  -0.745 0.456122    
date2024-04-16:treatmentPC1 -0.26135    0.26897  -0.972 0.331209    
date2024-01-12:treatmentPC2  0.37069    0.31011   1.195 0.231953    
date2024-01-23:treatmentPC2 -0.09502    0.27637  -0.344 0.730988    
date2024-02-07:treatmentPC2 -0.21252    0.29873  -0.711 0.476840    
date2024-02-21:treatmentPC2 -0.04440    0.29081  -0.153 0.878665    
date2024-03-06:treatmentPC2 -0.09098    0.30691  -0.296 0.766902    
date2024-03-22:treatmentPC2 -0.05387    0.29106  -0.185 0.853159    
date2024-04-04:treatmentPC2 -0.08550    0.29353  -0.291 0.770830    
date2024-04-16:treatmentPC2 -0.21145    0.30613  -0.691 0.489735    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)    -19.745   2462.968  -0.008    0.994
treatmentMC1    -4.190  19409.569   0.000    1.000
treatmentMC2    -4.202  19239.776   0.000    1.000
treatmentPC0    -4.212  19067.084   0.000    1.000
treatmentPC1    -4.184  19493.173   0.000    1.000
treatmentPC2    15.586   2462.969   0.006    0.995

simulationOutput <- simulateResiduals(fittedModel = fit_mg, n = 500)
plot(simulationOutput)

Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
: Fitting terminated with step failure - check results carefully

testUniformity(simulationOutput)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.0445, p-value = 0.4067
alternative hypothesis: two-sided

testOutliers(simulationOutput)

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 1, observations = 400, p-value = 1
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
 6.329252e-05 1.384977e-02
sample estimates:
frequency of outliers (expected: 0.00399201596806387 ) 
                                                0.0025 

testDispersion(simulationOutput)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.0908, p-value = 0.612
alternative hypothesis: two.sided

testZeroInflation(simulationOutput)

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 1.0941, p-value = 1
alternative hypothesis: two.sided

library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_mg, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") +
  theme_bw() +
  labs(
    x = "Date",
    y = "Mg",
    title = "Predicted Values of Mg",
    color = "Treatment"
  )

print(pred_effects)

# Predicted values of mg

treatment: MC0

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     48.98 | 33.19,  72.26
2024-01-12 |     59.67 | 41.83,  85.12
2024-01-23 |     67.20 | 48.77,  92.58
2024-02-07 |     55.67 | 39.96,  77.55
2024-02-21 |     60.49 | 43.37,  84.36
2024-03-06 |     68.08 | 48.83,  94.92
2024-03-22 |     79.82 | 57.83, 110.17
2024-04-04 |     85.52 | 61.87, 118.22
2024-04-16 |     74.36 | 53.75, 102.88

treatment: MC1

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     92.21 | 66.66, 127.56
2024-01-12 |     98.78 | 72.47, 134.65
2024-01-23 |     31.54 | 18.42,  54.01
2024-02-07 |     63.39 | 45.70,  87.93
2024-02-21 |     66.12 | 48.00,  91.10
2024-03-06 |     67.54 | 49.15,  92.81
2024-03-22 |     78.71 | 57.72, 107.34
2024-04-04 |     82.53 | 59.55, 114.38
2024-04-16 |     63.58 | 45.95,  87.99

treatment: MC2

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     62.99 | 44.42,  89.31
2024-01-12 |     85.50 | 62.61, 116.77
2024-01-23 |     95.53 | 70.42, 129.58
2024-02-07 |     73.32 | 53.58, 100.34
2024-02-21 |     86.61 | 63.86, 117.48
2024-03-06 |    100.09 | 74.39, 134.66
2024-03-22 |    117.84 | 88.03, 157.73
2024-04-04 |    116.16 | 86.53, 155.94
2024-04-16 |     70.77 | 50.38,  99.39

treatment: PC0

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     25.26 | 16.40,  38.90
2024-01-12 |     39.57 | 27.40,  57.16
2024-01-23 |     46.05 | 32.79,  64.68
2024-02-07 |     41.31 | 29.20,  58.45
2024-02-21 |     47.27 | 33.81,  66.09
2024-03-06 |     51.96 | 37.45,  72.11
2024-03-22 |     54.80 | 39.63,  75.76
2024-04-04 |     58.73 | 42.65,  80.86
2024-04-16 |     54.58 | 39.20,  75.99

treatment: PC1

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     44.01 | 28.78,  67.32
2024-01-12 |     52.44 | 36.28,  75.79
2024-01-23 |     49.82 | 34.42,  72.12
2024-02-07 |     40.04 | 28.05,  57.15
2024-02-21 |     45.29 | 32.09,  63.93
2024-03-06 |     53.25 | 38.18,  74.28
2024-03-22 |     61.05 | 44.13,  84.46
2024-04-04 |     63.28 | 45.15,  88.71
2024-04-16 |     51.46 | 35.63,  74.31

treatment: PC2

date       | Predicted |        95% CI
--------------------------------------
2023-12-27 |     34.13 | 20.75,  56.15
2024-01-12 |     60.24 | 41.87,  86.68
2024-01-23 |     42.59 | 29.66,  61.16
2024-02-07 |     31.37 | 21.44,  45.89
2024-02-21 |     40.33 | 28.26,  57.55
2024-03-06 |     43.32 | 29.98,  62.60
2024-03-22 |     52.71 | 37.58,  73.94
2024-04-04 |     54.72 | 38.69,  77.38
2024-04-16 |     41.95 | 29.13,  60.41

Adjusted for:
*     site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)

P

fit_p <- glmmTMB(
    p ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
    zi = ~ treatment,
    family = glmmTMB::lognormal(link = "log"),
    data = df_model
)

S(fit_p)

 Family: lognormal  ( log )
Formula:          
p ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation:     ~treatment
Data: df_model

      AIC       BIC    logLik -2*log(L)  df.resid 
   -807.8    -556.2     466.9    -933.8       338 

Random effects:

Conditional model:
 Groups            Name        Variance Std.Dev.
 site_block_id     (Intercept) 0.003353 0.05790 
 trt_site_block_id (Intercept) 0.007119 0.08438 
Number of obs: 401, groups:  site_block_id, 8; trt_site_block_id, 48

Dispersion parameter for lognormal family (): 0.0619 

Conditional model:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 -1.538415   0.151526 -10.153  < 2e-16 ***
date2024-01-12              -0.880465   0.265132  -3.321 0.000897 ***
date2024-01-23              -0.808236   0.196075  -4.122 3.75e-05 ***
date2024-02-07              -1.320144   0.235561  -5.604 2.09e-08 ***
date2024-02-21              -1.515269   0.306446  -4.945 7.63e-07 ***
date2024-03-06              -1.571143   0.262633  -5.982 2.20e-09 ***
date2024-03-22              -1.364637   0.268976  -5.073 3.91e-07 ***
date2024-04-04              -0.891213   0.251338  -3.546 0.000391 ***
date2024-04-16              -1.407929   0.265990  -5.293 1.20e-07 ***
treatmentMC1                -0.099689   0.195277  -0.511 0.609700    
treatmentMC2                -0.158669   0.190661  -0.832 0.405295    
treatmentPC0                -0.639774   0.262584  -2.436 0.014832 *  
treatmentPC1                 0.121799   0.199394   0.611 0.541303    
treatmentPC2                -0.009262   0.239372  -0.039 0.969134    
date2024-01-12:treatmentMC1  0.409739   0.339480   1.207 0.227447    
date2024-01-23:treatmentMC1 -0.091786   0.317385  -0.289 0.772433    
date2024-02-07:treatmentMC1 -0.051902   0.348987  -0.149 0.881774    
date2024-02-21:treatmentMC1  0.056781   0.392339   0.145 0.884929    
date2024-03-06:treatmentMC1  0.161877   0.373965   0.433 0.665111    
date2024-03-22:treatmentMC1  0.052483   0.418521   0.125 0.900205    
date2024-04-04:treatmentMC1 -0.772945   0.369608  -2.091 0.036505 *  
date2024-04-16:treatmentMC1 -0.220208   0.418549  -0.526 0.598804    
date2024-01-12:treatmentMC2  0.407385   0.314855   1.294 0.195707    
date2024-01-23:treatmentMC2 -0.268135   0.294249  -0.911 0.362162    
date2024-02-07:treatmentMC2 -0.014808   0.355481  -0.042 0.966772    
date2024-02-21:treatmentMC2  0.077202   0.399456   0.193 0.846749    
date2024-03-06:treatmentMC2  0.129368   0.406547   0.318 0.750325    
date2024-03-22:treatmentMC2 -0.216069   0.413588  -0.522 0.601374    
date2024-04-04:treatmentMC2 -0.234763   0.368105  -0.638 0.523630    
date2024-04-16:treatmentMC2  0.299460   0.371436   0.806 0.420116    
date2024-01-12:treatmentPC0  0.809856   0.386503   2.095 0.036140 *  
date2024-01-23:treatmentPC0  0.070908   0.350091   0.203 0.839494    
date2024-02-07:treatmentPC0  0.439496   0.394425   1.114 0.265163    
date2024-02-21:treatmentPC0  0.755476   0.442313   1.708 0.087634 .  
date2024-03-06:treatmentPC0  0.700833   0.419784   1.670 0.095016 .  
date2024-03-22:treatmentPC0  0.378590   0.428899   0.883 0.377397    
date2024-04-04:treatmentPC0  0.075187   0.388650   0.193 0.846601    
date2024-04-16:treatmentPC0  0.504392   0.389418   1.295 0.195235    
date2024-01-12:treatmentPC1  0.277727   0.324591   0.856 0.392207    
date2024-01-23:treatmentPC1 -0.516385   0.300661  -1.718 0.085887 .  
date2024-02-07:treatmentPC1 -0.113096   0.323280  -0.350 0.726459    
date2024-02-21:treatmentPC1 -0.029375   0.421446  -0.070 0.944432    
date2024-03-06:treatmentPC1 -0.095028   0.369877  -0.257 0.797243    
date2024-03-22:treatmentPC1 -0.285848   0.368584  -0.776 0.438027    
date2024-04-04:treatmentPC1 -0.542964   0.349546  -1.553 0.120341    
date2024-04-16:treatmentPC1 -0.212736   0.384762  -0.553 0.580330    
date2024-01-12:treatmentPC2  0.587871   0.356377   1.650 0.099030 .  
date2024-01-23:treatmentPC2 -0.017977   0.315887  -0.057 0.954618    
date2024-02-07:treatmentPC2  0.197177   0.348137   0.566 0.571137    
date2024-02-21:treatmentPC2  0.012480   0.418305   0.030 0.976199    
date2024-03-06:treatmentPC2  0.065331   0.383769   0.170 0.864826    
date2024-03-22:treatmentPC2  0.302133   0.397506   0.760 0.447212    
date2024-04-04:treatmentPC2  0.319921   0.349066   0.917 0.359403    
date2024-04-16:treatmentPC2  0.324149   0.383930   0.844 0.398506    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.53393    0.33244  -4.614 3.95e-06 ***
treatmentMC1  0.51228    0.43136   1.188    0.235    
treatmentMC2  0.63784    0.42535   1.500    0.134    
treatmentPC0  0.03843    0.45242   0.085    0.932    
treatmentPC1 -0.59769    0.51993  -1.150    0.250    
treatmentPC2 -0.42968    0.50305  -0.854    0.393    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

simulationOutput <- simulateResiduals(fittedModel = fit_p, n = 500)
plot(simulationOutput)

testUniformity(simulationOutput)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.087536, p-value = 0.004287
alternative hypothesis: two-sided

testOutliers(simulationOutput)

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 3, observations = 401, p-value = 0.2166
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
 0.001545489 0.021706793
sample estimates:
frequency of outliers (expected: 0.00399201596806387 ) 
                                           0.007481297 

testDispersion(simulationOutput)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 2.5953, p-value = 0.008
alternative hypothesis: two.sided

testZeroInflation(simulationOutput)

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 1.0005, p-value = 1
alternative hypothesis: two.sided

library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_p, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") + 
  theme_bw() +
  labs(
    x = "Date",
    y = "P",
    title = "Predicted Values of P",
    color = "Treatment"
  )